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Abstract 

A  general  approach  is  developed  for  the  rigorous  analysis  of  radiation,  scattering,  and  guid¬ 
ance  of  electromagnetic  fields  by  conducting  objects  of  arbitrary  shape  embedded  in  layered 
dielectric  media.  This  approach  is  based  on  the  mixed-potential  form  of  the  electric  field 
integral  equation,  which  is  amenable  to  the  existing,  well-established  numerical  solution 
procedures,  originally  developed  for  objects  in  homogeneous  space.  Numerical  results  are 
presented  for  surfaces  and  wires  that  penetrate  an  interface  between  dissimilar  media  and 
for  open  mi<  rostrip  transmission  lines  of  finite  thickness.  Computed  and  measured  data  are 
also  given  for  coax-fed  rectangular  and  triangular  microstrip  patch  antennas. 
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Chapter  1 
Introduction 


Simple  and  efficient  mclhou  of  muu^ats  (MOM)  [1]  procedures  have  recently  been  developed 
for  the  solution  of  the  electromagnetic  scattering  and  radiation  problems  involving  objects 
of  arbitrary  shape  [2,  3,  4,  5,  6,  7].  These  procedures  are  based  on  the  mixed-potential  form 
of  the  electric  field  integral  equation  (EFIE)  -  so  named,  because  it  involves  both  the  vector 
and  scalar  potentials,  the  former  expressed  in  terms  of  the  induced  current,  and  the  latter  in 
terms  of  the  induced  charge.  In  the  case  of  perfectly  electrically  conducting  (PEC)  objects, 
the  EFIE  is  more  general  than  the  magnetic  field  integral  equation  (MFIE)  [8],  since  it  is 
applicable  to  both  closed  and  open  surfaces  [9].  The  mixed-potential  EFIE  (MPIE,  for  short) 
is  preferable  to  several  other  possible  variants  of  the  EFIE,  because  it  only  requires  potential 
forms  of  the  Green’s  functions,  which  are  less  singular  than  their  derivatives  encountered 
in  other  forms  of  the  EFIE  [4].  In  the  MOM  technique  originally  developed  by  Rao  et 
al.  [5],  the  surface  of  the  PEC  object  is  modeled  in  terms  of  triangular  patches  and  specially 
designed  basis  functions  defined  on  pairs  of  adjacent  triangles  are  used,  which  yield  a  surface 
current  representation  free  of  line  or  point  charges  at  subdomain  boundaries.  This  technique 
also  employs  a  testing  scheme  in  which  the  derivatives  of  the  scalar  potential  are  in  effect 
replaced  by  finite  differences.  More  recently,  Schaubert  et  al.  [10]  extended  this  procedure  to 
volume  integral  equations  for  penetrable  bodies,  which  they  modeled  in  terms  of  tetrahedral 
elements. 

The  procedures  described  above  were  originally  developed  for  antennas  and  scatterers 
residing  in  a  homogeneous  space.  Aithough  this  restriction  is  not  severe  in  some  aerospace 
applications  where  the  effect  of  the  environment  can  be  neglected,  it  does  exclude  many 
problems  of  practical  interest  in  which  the  proximity  of  the  earth  must  be  taken  into  account. 
Indeed,  often  the  influence  of  the  ground  or  the  ocean,  which  in  many  cases  can  be  adequately 
represented  by  a  model  consisting  of  one  or  more  planar,  dielectric  layem,  is  the  dominant 
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effect  in  the  problem.  Microwave  and  millimeterwave  integrated  circuits  also  consist  of 
conductors  embedded  in  layered  dielectric  media. 

In  this  report  we  develop  MPIE  formulations  for  arbitrarily  shaped  objects  in  layered  me¬ 
dia  that  are  amenable  to  the  well-established  MOM  procedures,  such  as  those  implemented 
in  the  triangle-patch  code  [5,  7,  11]  and  its  thin-wire  counterpart,  MININEC  [12,  13].  Since 
a  great  amount  of  effort  went  into  the  development  of  these  codes  -  which  can  presently  only 
handle  objects  in  free  space  or  over  PEC  ground  -  it  is  desirable  to  have  a  formulation  that 
would  allow  one  to  easily  extend  them  to  the  case  of  bodies  in  stratified  media  Although, 
for  simplicity,  we  limit  attention  to  PEC  objects  and  surface  integral  equations,  our  formu¬ 
lation  can  easily  be  extended  to  dielectric  bodies  in  conjunction  with  either  volume  [10]  or 
surface  [6]  integral  equations  Unlike  in  free  space,  or  even  in  some  complex  environments 
with  PEC  boundaries,  it  is  not  a  trivial  task  to  formulate  an  MPIE  for  objects  of  arbitrary 
shape  in  layered  media,  because  in  such  environments  the  vector  and  scalar  potentials  are 
not  unique  [14]  and  the  scalar  potentials  of  point  charges  associated  with  horizontal  and 
vertical  dipoles  are  not,  in  general,  identical  [15].  Additional  complications  arise  when  the 
objects  penetrate  one  or  more  of  the  interfaces  between  dielectric  layers. 

The  advantages  of  the  MPIE  over  the  other  forms  of  the  EFIE  are  even  more  pronounced 
in  a.  layered  medium  than  in  free  space,  because  the  Green’s  functions  for  layered  media  com¬ 
prise  Sommerfeld-type  integrals  [16],  which  are  extremely  laborious  to  evaluate.  Since  the 
MPIE  involves  potential  forms  of  the  Green’s  functions,  the  Sommerfeld  integrals  it  requires 
converge  faster  than  those  present  in  any  other  form  of  the  EFIE  [4],  These  advantages  were 
previously  recognized  by  several  researchers.  Mosig  and  Gardiol  [17,  18]  have  applied  the 
MPIE  in  conjunction  with  the  Glisson  and  Wilton’s  [4]  MOM  procedure  to  a  rectangular 
microstrip  antenna.  Johnson  [19]  has  used  a  similar  approach  to  solve  the  problem  of  a  verti¬ 
cal  cylinder  penetrating  the  interface  between  contiguous  half-spaces.  Wilton  and  Singh  [20] 
have  applied  an  MPIE  to  analyze  a  periodic  array  of  slots  in  a  conducting  screen  backed 
with  a  layered  dielectric.  Michalski  et  al.  [21,  22]  have  used  an  MPIE  to  solve  the  problem 
of  a  horizontal  two-element  wire  array  above  and  on  opposite  sides  of  the  interface  between 
two  media.  As  was  pointed  out  in  [23,  24,  25],  the  success  of  these  efforts  can  be  attributed 
to  the  fact  that  the  structures  considered  could  only  support  either  vertical  or  horizontal 
components  of  the  current.  To  our  knowledge,  an  MPIE  for  arbitrarily  shaped  objects  in  a 
layered  medium  was  first  published  in  [24].  However,  it  was  assumed  in  that  paper  that  the 
antenna  or  scatterer  was  confined  to  a  single  layer.  In  a  two-dimensional  case,  an  MPIE  has 
been  derived  by  Xu  [26]. 


Numerous  papers  nave  been  published  on  the  subject  of  antennas  and  scatteiers  in  lay¬ 
ered  media,  but  -  with  the  exception  of  the  geophysics  literature  -  most  of  them  only  deal 
with  planar  geometries,  such  as  microstrip  antennas,  transmission  lines,  etc.  Since  the  em¬ 
phasis  here  is  on  objects  of  arbitrary  shape,  we  do  not  review  papers  in  this  cate'  ory,  to 
conserve  space.  The  problem  of  arbitrarily  shaped  thin  wires  which  are  near  to  or  penetrate 
an  interface  between  contiguous  half-spaces  has  been  solved  by  Burke  and  Miller  [27].  How¬ 
ever,  their  approach,  which  is  implemented  in  the  powerful  Numerical  Electromagnetics  Code 
(NEC)  [28],  is  not  easily  extendable  to  arbitrary  surfaces.  From  the  many  works  devoted  to 
electromagnetic  modeling  of  buried  inhomogeneities  in  the  context  of  geophysical  prospect¬ 
ing,  we  only  mention  the  recent  representative  papers  by  Holimann  [29]  and  by  Wannamaker 
et  al.  [30].  These  authors  use  the  volume  integral  equation  technique  in  conjunction  with  a 
rather  crude  -  but  entirely  adequate  in  the  quasi-static  regime  -  MOM  procedure  employing 
piece- wise  constant  current  expansion  and  point-matching  [1],  To  overcome  the  problems 
associated  with  the  singular  behavior  of  the  electric  Green’s  function,  Hohmann  [29]  em¬ 
ployed  a  mixed-potential  formulation,  but  only  to  the  primary  (or  whole-space)  component 
of  the  kernel;  the  part  comprising  the  Sommerfeld  integrals  was  left  in  the  slowly  convergent 
field  form.  Mention  should  also  be  made  of  the  work  by  Karlsson  and  Kristensson  [31], 
who  employed  the  extended  boundary  condition  integral  equation  in  conjunction  with  the 
T- matrix  approach  to  compute  the  field  scattered  by  obstacles  buried  in  a  stratified  ground. 
This  method,  however,  is  only  applicable  to  closed,  smooth  bodies  and  is  only  practical  for 
simple  shapes. 

The  outline  of  this  report  is  as  follows.  Chapter  2  contains  the  statement  of  the  problem 
and  a  general  discussion  of  various  forms  of  the  EFIE  in  a  layered  medium.  In  Chapter  3, 
alternative  forms  of  the  vector  potential  Green’s  function  are  derived  for  the  layered  medium. 
These  Green’s  functions  are  utilized  in  Chapter  4,  where  three  alternative  mixed-potential 
formulations,  referred  to  ac  jrmulations  A,  B  and  C,  are  developed  for  arbitrarily  shaped 
PEC  objects.  Formulation  C,  which  is  found  to  be  particularly  well  suited  for  the  application 
of  the  moment  method,  is  specialized  in  Appendix  A  to  two  particularly  important  situations, 
where  the  object  resides  in  contiguous  half-spaces  or  is  partially  embedded  in  a  grounded 
dielectric  slab.  In  Chapter  5,  the  MOM  procedures  introduced  in  [4,  5,  7]  are  adapted  to 
solve  the  MPIE  based  on  Formulation  C  of  Chapter  4.  For  surfaces  of  arbitrary  shape, 
the  triangle-patch  model  [5]  is  employed,  while  thin-wire  structures  are  modeled  by  straight 
tubular  segments  [7].  The  wire-to-patch  junction  model  developed  by  Hwu  et  al.  [7]  has 
been  adopted  for  the  analysis  of  coax -fed  microstrip  patch  antennas.  Formulation  C  is  also 


specialized  in  Chapter  5  to  the  case  of  an  open  transmission  line  consisting  of  an  infinite 
conductor  of  arbitrary  cross-section  partially  embedded  in  a  grounded  dielectric  slab.  In 
the  numerical  procedure,  the  conductor’s  cross-section  profile  is  modeled  by  straight  line 
segments  [4].  In  Chapter  6,  the  numerical  procedures  used  to  evaluate  the  Sommerfeld 
integrals  are  discussed.  In  Chapter  7,  sample  numerical  results  are  presented  for  several  cases 
of  interest  involving  PEC  plates,  cylinders,  thin-wire  antennas,  coax-fed  microstrip  patch 
antennas,  and  open  microstrip  lines.  Where  possible,  comparisons  are  made  with  measured 
data  or  with  computed  results  obtained  by  other  researchers.  Finally,  in  Chapter  8,  we  draw 
conclusions  and  make  recommendations  for  future  work. 
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Chapter  2 
P  reliminar  ies 


2.1  Statement  of  the  Problem 

Consider  a  medium  consisting  of  n  +  1  dielectric  layers  separated  by  n  planar  interfaces 
parallel  to  the  xj/-plane  of  a  Cartesian  coordinate  system  and  located  at  z  =  z\,  l  =  1,2,. .  .,n, 
as  illustrated  in  Fig.  2.1.  The  medium  of  the  z th  layer  is  characterized  by  permeability  pi 
and  permittivity  e^,  which  may  be  complex  if  the  medium  is  lossy.  Let  the  PEC  object  (or 
collection  of  objects)  in  Fig.  2.1  occupy  p  layers  with  indices  /j,  l2, . . . ,  /p,  where  1  <  p  <  n-fl. 
For  later  convenience,  define  the  ordered  set  of  indices  L  —  {lx,  l2, . . . ,  /p},  in  which  <  h t+i- 
Let  Si  denote  the  surface  of  the  object(s)  in  the  zth  layer  and  let  h,  be  a  unit  vector  normal 
to  Si.  The  quantity  of  interest  is  the  surface  current  density  J(r)  excited  on  the  object(s) 
by  a  given  time-harmonic  incident  electric  field  E,nc.  The  time  variation  is  assumed  and 
suppressed. 


2.2  Electric  Field  Integral  Equation  (EFIE) 

The  EFIE  for  the  current  density  J  on  the  surface  S  of  the  PEC  object(s)  embedded  in  a 
layered  medium  is  obtained  by  enforcing  the  boundary  condition  [32] 


-hm  x  E*m(r)  =  nm  x  E'”c(r),  r  on  Sm,  me  L 


(2.1) 


where  r  is  the  position  vector  defined  with  respect  to  a  global  coordinate  origin,  E™c  is  the 
“incident”  electric  field  (i.e.,  the  field  in  the  absence  of  the  object)  in  the  mth  layer,  and  Eam 
is  the  scattered  field  in  the  mth  layer.  For  the  structure  of  Fig.  2.1,  E*m  and  the  scattered 
magnetic  field  H3m  can  be  expressed  as 


-E'Jr)  =  £  +  V<T'(r)] 


(2.2) 


x€L 
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(2.3) 


W=-VxE^‘(r) 

m  i€L 

where  Am‘  is  the  magnetic  vector  potential  in  the  mth  layer  due  to  the  current  density  J  in 
the  ith  layer,  and  is  given  as 


>)  =  jG7(r\r')-J(r')dS' 


and  4>m,(r)  is  the  corresponding  scalar  potential  which  is  related  to  Ami(r)  through  the 
Lorentz  condition 

<T‘(*-)  =  ^V-Am‘(r)  (2.5) 

m 

where  k ^  =  u>2em/im.  In  (2.4),  is  the  dyadic  Green’s  function  which  represents  the 
magnetic  vector  potential  in  region  m  due  to  a  unit-strength,  arbitrarily  oriented  current 
dipole  in  region  i.  G™1  can  be  found  by  solving  the  inhomogeneous  Helmholtz  equation 


(V2  +  kl)  I  r')  =  ~Hm  L6  (**  —  O  (2.6) 

where  £  is  the  idemfactor,  subject  to  the  condition  that  the  tangential  component  of  Eam 
and  H3m  be  continuous  across  the  interfaces  between  dielectric  layers.  As  is  well  known  [16], 
for  a  horizontal,  say,  x-directed  dipole,  two  components  of  the  vector  potential  are  required 
to  satisfy  the  boundary  conditions  at  the  interfaces.  Traditionally  [33],  the  z  component  has 
been  selected  in  addition  to  the  x  component.  The  Green’s  function  in  this  case  takes  the 
form  [24] 

=  (x  x  +  y  y)  G™  +  »  x  G™  +  z  y  G™  +  z  z  G™'.  (2.7) 

However,  one  may  as  well  postulate  the  y  component  of  the  vector  potential  to  accompany 
the  primary  x  component  [14].  This  strategy  leads  to  an  alternative  form  of  the  Green’s 
function, 

fi?  =  *  *  Gt  +  yyG™  +  (xy  +  yx)  G™  +  z  zG™.  (2.8) 

In  (2.7)  and  (2.8)  G™  denotes  the  x-component  of  Am'  due  to  an  x-directed,  unit-strength 
electric  dipole,  G™  the  x-component  of  .A”1'  due  to  a  similar,  y-d irected  dipole,  etc.  We 
note  that,  except  for  G£',  the  corresponding  components  in  (2.7)  and  (2.8)  are  different,  even 
though  the  same  notation  is  used.  The  entries  of  the  dyadics  in  (2.7)  and  (2.8)  are  developed 
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in  the  next  chapter,  and  are  shown  to  comprise  improper  integrals  of  the  Sommerfeld  [16] 
type. 

Still  other  Green’s  functions,  in  addition  to  those  in  (2.7)  and  (2.8)  are  possible.  For 
example,  one  can  obtain  another  form  of  by  postulating  that  the  y  and  2  components  of 
the  vector  potential,  instead  of  the  x  and  2  or  x  and  y  components,  accompany  an  x-directed 
dipole.  However,  this  and  other  forms  of  Q.™  are  less  attractive  for  our  purpose  than  (2.7) 
and  (2.8),  and  are  not  considered  here. 

Substitution  of  (2.2)  into  (2.1)  yields 

nm  x  £  [juAmi(r)  +  V<T‘(r)]  =  hm  x  E'™(r),  r  on  Sm,  rn  €  L  (2.9) 

x€L 

which,  by  invoking  (2.4)  and  (2.5),  can  be  further  transformed  to 
x  (VV  ■  +Ai)  £  /£”’(,  | ,')  ■  J(r')  iff 

771  I  €  Lt  g . 

=  nm  x  E'*c(r),  r  on  Sm,  m  G  L.  (2.10) 

This  equation  is  referred  to  as  vector-potential  EFIE  [34],  since  it  only  involves  the  magnetic 
vector  potential.  Although  (2.10)  has  often  been  used  in  MOM  analyses  of  straight  wire 
antennas  or  planar  scatterers,  the  presence  of  the  mixed  tangential  derivatives  makes  it  less 
suitable  for  objects  of  arbitrary  shape.  By  introducing  the  differential  operator  under  the 
integral  sign  in  (2.10),  we  may  obtain  another  form  of  the  EFIE,  in  which  the  dyadic  kernel 
is  the  Green’s  function  for  the  electric  field.  However,  this  EFIE  is  not  attractive  because  the 
kernel  is  highly  singular,  which  makes  the  evaluation  of  the  integrals  required  by  the  MOM 
procedure  difficult  when  the  observation  point  is  within  the  integration  interval  [35].  Also,  the 
required  differentiation  of  the  Sommerfeld- type  integrals  adversely  affects  their  convergence. 
These  difficulties  can  be  avoided  if  only  one  of  the  operators  nabla  is  introduced  inside  the 
integral  in  (2.10)  and  then  transferred,  by  a  series  of  transformations,  to  act  on  the  current. 
The  result  is  the  mixed-potential  EFIE  discussed  below. 

2.3  Mixed-Potential  EFIE  (MPIE) 

We  note  that  (2.9)  would  be  in  the  desired  mixed-potential  form  if  the  scalar  potential 
were  expressed  in  terms  of  the  surface  charge  density  q{r).  With  this  goal  in  mind,  we 
substitute  (2.4)  into  (2.5)  and  introduce  the  operator  nabla  under  the  integral  sign  (this  step 
can  be  justified  [36,  37])  to  obtain 
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m  S, 


(2.11) 


Obviously,  our  objective  would  be  achieved  if  we  transferred  the  divergence  operator  to  act 
on  the  current,  in  view  of  the  equation  of  continuity,  V  J  =  —juq.  It  is  shown  below  that 
this  can  only  be  accomplished  if  a  scalar  function  G™  can  be  found,  such  that 

.  gj‘(r  I  r')  =  J-- V'G"(t  |  r').  (2.12) 

In  a  homogeneous  medium,  where  GJ'1  may  be  interpreted  as  the  Green’s  function  for  the 
scalar  potential,  this  is  a  quite  trivial  task.  If  the  medium  is  stratified,  however,  G sat¬ 
isfying  (2.12)  does  not  in  general  exist  [24],  which  can  be  attributed  to  the  fact  that  the 
scalar  potentials  of  point  charges  associated  with  horizontal  and  vertical  current  dipoles  in 
a  layered  medium  are  in  general  different  [15].  Hence,  in  order  to  achieve  our  goal,  we  follow 
the  procedure  proposed  in  [23,  24],  and  introduce  a  scalar  function  K™  and  a  vector  function 
Pmx  according  to 

pr  V  ■  Gf(T  I  r')  =  ^V'A7'(r  I  r')  +  jwPm'(r  |  r').  (2.13) 


One  should  note  that  in  the  above  the  choice  of  K™  and  Pmt  is  not  unique. 

Upon  substituting  (2.13)  into  (2.11)  and  using  a  vector  identity  (p.  487,  [38])  and  the 
Gauss  theorem  [38,  p.  503],  we  can  express  the  scalar  potential  as 


r\r)  =  J  I<™(r  I  r')  q(r')  dS'  +  ju>  J  Pmi(r  |  r')  •  J(r')  dS' 

S,  Si 

\r')J(r')  ■  itidC1  —  j  K?(r  \  r')  J(r')  •  dC' 


1 

+  1 — 


(2.14) 


where  Ci  and  Ci-\  are  the  contours  formed  by  the  intersection  of  the  surface  5,-  with  the  inter¬ 
faces  at  z  =  Zi  and  z  =  Zi-\,  respectively,  and  it,  and  u,_i  are  the  unit  vectors  perpendicular 
at  r'  to  Ci  and  Ci- 1,  respectively,  in  the  planes  tangent  to  S;  (Fig.  2.2).  Substituting  (2.14) 
into  (2.9),  we  finally  obtain 


Hm 


jut  I  K™'(r  |  r')  •  J(r')  dS'  +  V 
5, 

fK?(T\r')J(r')-u,dC'-  j 

Ci  i 


/  Kf(r\r')q(r')dS' 

s, 

K?i(r\r')J(r')-ui-ldC' 


=  nm  x  E'™(r),  r  on  5m,  m  g  L 


(2.15) 
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(M+t,  *i* i 


'-'4+1 
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where  we  have  introduced  the  dyadic  kernel 


|  r')  =  G™(r  \  r')  +  V  Pmi{r  |  r').  (2.16) 

We  note  that  (2.15)  would  be  in  the  desired  mixed-potential  form  [4,  5]  if  it  were  not  for 
the  presence  of  the  term  contributed  by  the  contour  integrals,  which  occur  when  the  object 
penetrates  one  or  more  of  the  interfaces.  In  Chapter  4  we  show  that  with  a  proper  choice 
of  G1?1  and  K ™  in  (2.13)  the  contour  integrals  cancel  out.  We  note,  however,  that  even  if  a 
formulation  is  chosen  in  which  the  contour  integrals  persist,  the  MOM  procedures  developed 
in  [4,  5]  can  be  extended  to  accommodate  these  terms.  As  a  matter  of  fact,  the  “correction 
term”  V  Pm>  could  also  be  handled  in  this  manner,  instead  of  being  incorporated  into  the 
vector  potential  kernel  via  (2.16). 


Chapter  3 


Derivation  of  Green’s  Functions  for 
the  Magnetic  Vector  Potential 


The  dyadic  Green’s  function  G™*(r  |  r')  in  (2.4)  is  the  vector  potential  in  region  m  due  to 

•  a  unit-strength,  arbitrarily-oriented  Hertzian  dipole  in  the  ith  layer.  The  direct  solution 
of  (2.6)  for  G™‘  is  extremely  tedious,  since  various  components  of  this  dyadic  couple  in  the 
boundary  conditions  at  the  interfaces  between  dielectric  layers.  Therefore,  in  this  chapter,  we 
derive  the  Green’s  function  by  means  of  the  Fourier  transformation,  which  in  effect  reduces 

$  the  original  problem  to  that  of  solving  an  equivalent  transmission-line  network  along  the  2 

coordinate. 

3.1  Reduction  of  Maxwell’s  Equations  to  Transmis- 

•  sion  Line  Equations  for  Sources  in  Layered  Media 

For  the  geometry  of  Fig.  3.1a,  we  are  interested  in  the  electromagnetic  field  ( E,H )  every¬ 
where  in  the  medium  due  to  prescribed  distributions  of  time-harmonic  electric  and  magnetic 

•  currents  J  and  M ,  respectively.  Since  the  structure  is  infinite  along  the  x  and  y  coordinates, 
we  can  simplify  the  analysis  by  utilizing  the  “shifted”  Fourier  transform  pair 

T{f{x  -  x',y  -  y')}  =  f(kx,  ky) 

oo  oo 

•  =  /  J  f{x-x',y -y'W^-^+^-^Uxdy  (3.1) 

—  CO  — oo 

F~l{J{kx,ky)}  =  f(x  -x',y-  y') 

OO  OO 

•  =Wf!  j  (3.2) 

—  oo  — oo 
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ZfH  22 


(b) 


Figure  3.1:  (a)  Plane-stratified  dielectric  medium  and  (b)  its  transmission-line  network  rep¬ 
resentation. 
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where  ( x',y ')  locates  the  projection  of  the  source  point  on  the  xy- plane  (Fig.  3.2).  The 
inverse  Fourier  integrals  (3.2)  can  be  conveniently  expressed  in  terms  of  the  Sommerfeld 
integrals  [16]  of  the  form 


(3.3) 


by  changing  to  polar  coordinates  in  both  the  transform  and  space  domains  according  to 


x  —  x  =  £  cos  £,  y  —  y  =  £  sin  (  (3.4) 

kx  =  kpcosa,  ky  =  kpsina  (3.5) 


where  (Fig.  3.2) 


(  =  \/(x  -  *')2  +  {y  -  y')2<  C 


kp  =  \Jk2x+  a  —  arctan 


=  arctan 


ky 

K 


(3.6) 


(3.7) 


and  where  Jn  is  the  Bessel  function  of  order  n.  Hence,  using  (3.3).  we  can  express  (3.2)  as 


x-‘  {/(*,)}  =  4- So  [/(*,)  . 

(3.8) 

We  also  find  that 

x-'  {;*,/(*,)}  =  [/(*.)] 

(3.9) 

X-'UKHK))  =  s-~S,  [/(*.)! 

(3.10) 

=  -T  {cos2CS2[/(*,)]  -  Sb[*.7(M]} 

(3.11) 

=  ~{cos2C52[/(  *,)]  +  Sb[*.7<*,)l} 

(3.12) 

T~'  {kj;vs{kp)}  =  ~4-  sin2(  52[</(^)]. 

(3.13) 
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xy  -  plane 


projection  of  the 
observation  point 


Figure  3.2:  Definition  of  the  distance  £  and  angle  (. 
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These  relationships  will  be  used  in  the  formulation  of  vector  potential  Green's  functions.  We 
list  them  here  for  later  convenience. 

Since  the  structure  of  Fig.  3.1a  can  be  regarded  as  a  uniform  cylindrical  waveguide  (of 
infinite  cross-section)  with  z  as  the  longitudinal  (axial)  coordinate,  it  is  fruitful  to  decompose 
the  field  and  source  vectors  into  their  transverse  and  longitudinal  parts,  be.,  we  can  write  E  = 
Et  +  zEz,  H  —  H,  +  zHz ,  etc.  The  solution  of  the  field  equations  for  E  and  H  is  facilitated 
by  eliminating  the  dependent  longitudinal  components  from  the  Maxwell’s  equations  [39.  40]. 
Hence,  by  eliminating  Ez  and  H ,  and  invoking  the  Fourier  transformation  (3.1),  we  arrive 
at 


~~Et=j^  U 


d  fr  _  (r  VtVt 


+  —  tj^t  j  ■  (Ht  X  z)  +  Mte  X  z 
)-(*x£()  + 


Z  X  J 


(3.14) 

(3.15) 


where  the  equivalent  transverse  electric  and  magnetic  current  distributions  are  given,  respec¬ 
tively.  as 


J te  —  J l  Z  X 


and 


TJ  ^  .  VtJz 

Mte  =  Mt  +  z  x  — - . 

jut 


(3.16) 


(3.17) 


X  is  a  unit  dyadic,  and  V(  =  —jkxx  —jkyy  is  the  transform  domain  counterpart  of  the  trans¬ 
verse  operator  nabla.  The  longitudinal  field  components  are  derivable  from  the  transverse 
components  as  [40] 


• 

e2  =  —  {v(  . 

(Ht  x  z)  —  1} 

(3.18) 

]UJt  *• 

H,  =  —  {V( 

•  (z  x  Et)  -  Mz\ 

(3.19) 

\  /  ) 

• 

Equations  (3.14)-(3. 19)  hold  for  each  layer  of  the  stratified  medium  (Fig. 
subscript  has  been  suppressed  in  these  equations  for  simplicity. 

3.1a).  The  layer 

Since  any  vector  can  be  represented  as  a  sum  of  two  parts,  one  of  which 

is  solenoidal  and 

the  other  irrotational, 

we  may  express  the  vectors  Et  and  Ht  x  z  as  [40] 

• 

II 

1 

<\ 

<■» 

l 

■  VtVh  x  z 

(3.20) 
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Htxz  =  -Vtr~VtIh  X  z . 


(3.2i : 


Using  these  equations  in  (3.18)  and  (3.19),  we  can  express  the  longitudinal  field  components 


as 


e,  =  —  (V  r  - 1) 

ju>e  v  ' 

Hz  =  - —  ( kp2Vk  -  A#,) 

11. Ml  \  ' 


(3.22) 

(3.23) 


where  kp 2  =  k2  +  Equations  (3.20)— (3.23 )  indicate  that  in  source-free  regions,  ( Ve ,  Ie) 
and  (Vh,  Ih)  generate,  respectively,  field  transverse  magnetic  (TM)  and  transverse  electric 
(TE)  to  z.  To  determine  the  scalar  functions  Ve,  le,  Vh,  and  Ih ,  we  use  (3.20)  and  (3.21) 
to  eliminate  Et  and  Ht  from  (3.14)  and  (3.15).  As  a  result,  we  obtain  the  equations  [40] 

dVq 


dz 

ill 

dz 


=  jkzZq  Iq  +  vq 
=  jkzYq  Vq  +  iq 


(3.24) 

(3.25) 


where  the  superscript  “q”  stands  for  either  “e”  or  “A”  and  where  k]  =  k 2  —  kp2.  The  modal 
impedances  Zq  (admittances  Yq)  are  defined  as 


7e  =  —  =  — 
Ye  use 


Zh  =  _L  _  ^ 


V'1  A:, 

and  the  source  functions  vq  and  iq  are  given  by  [40] 

u'  =  -— +  ~(zx 

jue  k2  v  > 

•"  = 


(3.26) 

(3.27) 


Kn 


Mz  1  /.  ~  \  x 

- - rr  (2  x  v0  •  Jt- 

kp  v  > 


(3.28) 

(3.29) 

(3.30) 

(3.31) 
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Equations  (3.24)  and  (3.25)  have  the  form  of  transmission-line  equations  for  the  voltage 
Vq  and  current  Iq.  To  the  /th  layer  of  the  stratified  medium  (Fig.  3.1a)  there  corresponds 
a  transmission  line  section  with  characteristic  impedance  Zf  and  propagation  constant  kzr 
Hence,  the  layered  medium  of  Fig.  3.1a  may  be  represented  by  two  transmission-line  net¬ 
works,  which  may  be  referred  to  as  the  TM  and  TE  networks,  each  of  the  form  shown  in 
Fig.  3.16.  The  TM  network  has  characteristic  impedances  Zf  and  sources  (ve,  ie),  while 
the  TE  network  has  characteristic  impedances  Zj1  and  is  driven  by  ( vh ,  ik).  Depending 
on  the  sources  (cf.  (3.28)— (3.31)),  the  TM,  TE,  or  both  networks  may  be  excited.  These 
network  problems  must  be  solved  subject  to  the  condition  that  the  voltages  and  currents 
be  continuous  across  the  interfaces  between  contiguous  line  sections.  These  conditions  are 
the  consequence  of  the  continuity  of  the  transverse  components  of  the  electric  and  magnetic 
fields  across  the  interfaces  between  adjacent  dielectric  layers  (Fig.  3.1a). 

The  solution  of  the  transmission-line  equations  (3.24)  and  (3.25)  is  facilitated  if  one 
appeals  to  the  principle  of  superposition  and  considers  the  effect  of  the  voltage  and  current 
sources  separately.  Hence,  with  only  voltage  sources  present  ( iq  =  0),  (3.24)  and  (3.25) 
reduce  to 

Iq  =  jkzYqvq  (3.32) 


Vq  =  -  ——p 
j  kz  dz  ' 


(3.33) 


From  these  equations,  the  impedances  seen  looking  to  the  left  and  to  the  right  at  a  point  z 
on  the  transmission  line  can  be  expressed,  respectively,  as 

J,^  =  ~Wi  =  W,TPP(z)  (3'34) 

and 

^-m=-w,Tpp^  (3-35) 

Similarly,  with  only  current  sources  present  ( vq  =  0),  (3.24)  and  (3.25)  reduce  to 


l/?  =  jkzZqiq 


(3.36) 


Yq  d 

/<? 

j  kz  dz 


(3.37) 
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The  admittances  seen  looking  to  the  left  and  to  the  right  at  a  point  z  on  the  transmission 
line  can  be  in  this  case  expressed,  respectively,  as 


and 


<—  I«(z)  Yq  d 


I<(z)  Yq  d 


(3.38) 


(3.39) 


Once  the  transmission-line  equations  (3.32),  (3.33),  (3.36)  and  (3.37)  are  solved  for 
(Ve,/e)  and  ( Vh,Ih ),  we  determine  ( E,H )  from  (3.20)-(3.23).  Finally,  the  inverse  Fourier 
transformation  (3.2)  is  employed  to  recover  ( E,H ).  The  solution  of  (3.32),  (3.33),  (3.36) 
and  (3.37)  for  the  stratified  medium  of  Fig.  3.1a  is  discussed  in  the  next  section. 


3.2  Solution  of  Transmission-Line  Equations 

It  was  shown  in  the  last  section  that  the  layered  medium  of  Fig.  3.1a  can  be  represented  by 
the  transmission-line  network  of  Fig.  3.16.  The  voltage  and  the  current  in  each  transmission¬ 
line  section  are  governed  by  (3.24)  and  (3.25),  or  equivalently,  by  (3.32),  (3.33),  (3.36), 
and  (3.37).  In  solving  these  equations,  it  is  helpful  to  define  for  each  interface  the  impedances 
and  voltage  reflection  coefficients,  as  illustrated  in  Fig.  3.3.  (For  notational  simplicity,  we 
suppress  the  superscripts  “<?”.)  The  solution  of  (3.32),  (3.33),  (3.36),  and  (3.37)  is  simplified 
by  first  considering  point-source  excitations  and  then  obtaining  the  total  response  by  use 
of  the  superposition  theorem  [40].  Without  loss  of  generality,  we  assume  that  the  sources 
are  located  in  the  ith  section  of  the  transmission  line  (which  corresponds  to  the  ith  layer  in 
Fig.  3.1a)  and  that  the  impedances  Z,(zi)  =  Z,  and  Z,(z;_i)  =  Z,_!  (or  their  inverses  Yt 
and  T,_i)  are  specified.  Hence,  if  we  denote  the  current  /  in  the  ith  line  section  due  to  a 
point  voltage  v  =  —  6(z  —  z ')  in  the  same  section  as  Gjt,  we  obtain  from  (3.32) 


(3.40) 


This  equation  must  be  solved  subject  to  the  boundary  conditions  (cf.  (3.34)  and  (3.35)) 
Z,  d 


Z,  = 


jkZi  dz 


In  Gl(z,z'),  z  —  Z{ 


Z,-i  =  —■~-^r-\nGIll(z,z),  z  =  Z,_!. 
j  dz 


(3-41) 


(3.42) 
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i 

(b) 


Figure  3.4:  (a)  Voltage  source  and  (b)  current  source  in  the  zth  transmission-line  section. 
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Equations  (3.40)-(3.42)  can  be  represented  by  the  network  shown  in  Fig.  3.4a. 

In  a  like  manner,  if  we  denote  the  voltage  V  in  the  zth  line  section  due  to  the  point 
current  i  =  —6(z  —  z ')  in  the  same  section  as  Gl,  we  obtain  from  (3.36) 


£_ 

dz2 


Gl{z,z,)  =  -jk'iZi6{z-z'). 


(3.43) 


The  boundary  conditions  in  this  case  are  (c/.  (3.38)  and  (3.39)) 

Yi  =  -jrj-^GUz’  z'),  Z  =  Zi 
jk2idz 

Y i~i  =  — ^-j-ln  ■(*,«'),  z  =  z,-u 
J  kzi  GZ 


(3.44) 

(3.45) 


Equations  (3.43)-(3.45)  have  the  network  representation  shown  in  Fig.  3.46. 

We  first  turn  attention  to  (3.43);  the  solution  of  (3.40)  will  follow  by  duality.  Hence, 
solving  (3.43)  subject  to  the  boundary  conditions  in  (3.44)  and  (3.45),  we  obtain 


,~jkz,\z-z'  I 


+  Q*W) 


(3.46) 


with 


QY(z,z')  =  ~  +  T",_ie-jfc»[22-‘-(2+z')] 

Is  i  ^ 

+2Tl%_1e-W'  cos  [*„•(*“*') 


(3.47) 


A  =  1  -  r,  ri_ie"i2¥’<  (3.48) 

where  A  =  kzid{,  d<  =  z,_ i  —  z,-.  The  reflection  coefficients  T  *  and  T  *  (see  Fig.  3.3)  are 
given  by,  respectively, 


I '  Z  k  Z k+ 1  to 

*■  k  —  ^  j  ’  5  ^ 

Zk  +  Zk+ 1 


and 


V  Zk-  Zk  ,  ,  n 

*  k  —  ^ k  —  1  ?  2,  •  •  •  ,  71 

Zk  +  Zk 


(3.49) 


(3.50) 


where  r0  =  0  =  rn+i. 

The  yet  unspecified  impedances  Zk  and  Z k  can  be  obtained  from  the  analysis  of  the 
transmission-line  network  of  Fig.  3.16.  It  can  be  shown  that  these  impedances  are  given  by 
the  recursion  formulas 


and 


rp  rp  ^  k—l  J  ^  k  ^  k  i  0  9 

^ k  —  " k  ~  ^  k  * —  2, 3, * •  • ,  n 

Zk  +  j  z  k-\tk 


(3.51) 


V  7  Z*+1  +  i^+l^+l  ,  1  O  1 

"  k  ^k+1  , _  j  ^  ^  1?TZ  2,  *  ‘  ,  1 


(3.52) 


■2^/t+i  +  j  Z k+ih+i 

where  Z\  =  Zu  Zn  =  Zn+ 1,  and  4  =  tan  V»jt- 

From  (3.46),  one  can  compute  the  voltage  anywhere  in  the  ith  line  section,  including  the 
terminals,  due  to  a  unit  current  source  in  the  same  section  (cf.  Fig.  3.46).  The  voltage  in 
any  other,  say  mth  section  due  to  this  source  can  be  easily  determined  from  the  analysis  of 
the  transmission-line  network  of  Fig.  3.16  as 

GYiizi-uz'yf^z),  i  -  1  >  m  >  1 

Gvmt(z,z')=l  ^  (3.53) 

GYi(zi,  z')*T%t(z),  n  +  l>m>i  +  l 

where  the  voltage  transfer  functions  T^t  and  T ^  are  given,  respectively,  as 

P-jkZm(z-zm)  (  .  ^  i-2  fl  +  r c“^*+» 


= 


1  + 


1  +  Ym-Xe-W” 


n 


and 


TlA*)  = e 


~jkZm  {zm  —  l  ~~z) 


J  _j_  p  g*“i2fcxm(2— 2m) 


k=m  1  -f  r  +  i 

tn— 1  |  1  +  ‘f  b  £-J*» 


-(3.54) 


}•  n  - 

1  4=i+l  1 


+  f 


(3.55) 


1  +  Tm^‘2'Prn 

This  completes  the  solution  of  (3.43).  Equation  (3.40)  can  be  solved  by  a  dual  procedure. 
However,  this  is  not  necessary:  we  can  obtain  G G  from  G by  replacing  in  the  latter  all 
impedances  Z,  by  their  reciprocals  Yj  (as  a  result  of  this  operation,  all  reflection  coefficients 
change  signs,  as  is  evident  from  (3.49)-(3.50)).  The  functions  G!m%  and  G^  will  be  referred 
to,  respectively,  as  the  current  and  voltage  transmission-line  Green's  functions. 


3.3  Derivation  of  Dyadic  Green’s  Functions  for  the 
Magnetic  Vector  Potential  in  a  Plane- Stratified 
Medium 

The  magnetic  vector  potential  Amx  in  the  mth  layer  of  the  stratified  medium  (Fig.  3.1a)  due 
to  a  current  distribution  J  occupying  a  volume  V [  in  the  ith  layer  can  be  obtained  from 


Am\v)  =  jg™(r\r')-  J(r')dV' 


(3.56) 


where  G is  the  dyadic  Green’s  function  for  the  magnetic  vector  potential,  which  can  be 
expressed  either  in  the  traditional  form  (2.7),  or  the  alternative  form  (2.8).  To  determine 
the  elements  of  the  dyadic  in  (3.56),  we  appeal  to  the  transmission-line  analogy  and  make 
use  of  the  transmission-line  Green’s  functions  developed  in  Sections  3.1  and  3.2.  Hence, 
from  (3.29)  and  (3.31)  we  determine  that  an  x-directed  dipole  with  J  =  J  t  =  x6(z  —  z') 
gives  rise  to  transmission-line  sources  ie  and  ih  given,  respectively,  as 


<■(*)  -  -  *0 

Kp 


(3.57) 


H*)  =  -Jrr<5(z  -  z') 

Kp 


(3.58) 


where  we  assume  that  the  sources  are  in  the  ith  layer,  t.e.,  2,  <  z'  <  2,_!.  Obviously,  both 
the  TM  and  TE  transmission-line  networks  are  excited  in  this  case.  The  current  generators 
given  by  (3.57)  and  (3.58)  excite,  respectively,  the  voltages  and  K^i  in  the  mth  sections 
of  the  transmission-line  network.  From  Section  3.2,  we  find 


KM  =  rr<%(*,*') 


(3.59) 


KM  =  t?G%(z,z'). 

Kp 

For  a  ^-directed  dipole  with  J  =  Jt  =  y&{z  —  z'),  (3.29)  and  (3.31)  yield 


(3.60) 


ie(z)  = -J-ri6(z  -  z') 

Kp 


(3.61) 


lh(z)  =  +TT(5(2  -  2')- 


(3.62) 


These  sources  excite  the  voltages  (cf.  Section  3.2) 
KM=TiG%(z,z') 


(3.63) 
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ViW  =  -^G5s(*,*')  (3.64) 

ftp 

respectively. 

Finally,  for  a  vertical  dipole  with  J  =  zJ2  =  z6(z  —  z '),  we  obtain  from  (3.28) 

ue(2)  =  —  - — 6(z  —  z).  (3.65) 

JUti 

Hence,  only  the  TM  transmission-line  network  is  excited  in  this  case.  The  current  in  the 
mth  section  due  to  the  voltage  source  give  in  (3.65)  is  easily  found  as  (cf.  Section  3.2) 


]uex 


(3.66) 


We  now  must  relate  these  transmission-line  voltages  and  currents  to  the  components  of 
the  magnetic  vector  potential.  First,  we  turn  our  attention  to  a  vertical  dipole  source,  which 
gives  rise  to  only  the  2-component  of  the  vector  potential.  Outside  the  source  (2  ^  2'),  the 
2-component  of  the  electric  field  can  be  expressed  in  terms  of  A™1  as 


Fm'  —  ^  4-  k 2  )Ami  —  2  Ami 

^ *  ~  £2  +  km)A*  -  ■ 

On  the  other  hand,  from  (3.22)  we  have 


(3.67) 


E™  = 


-k  2r 

KP  1m‘ 


(3.68) 


Eliminating  E ™l  between  (3.67)  and  (3.68),  we  obtain 
Ami  —  11  Te 


(3.69) 


It  can  be  shown  that  this  solution  also  holds  in  the  source  region  (2  =  z').  If  the  current 
in  (3.69)  is  that  generated  by  a  2-directed  dipole,  then  A™1  =  G™.  Hence,  substituting 
from  (3.66)  into  (3.69),  we  obtain 


1  k-m  p^ir 


/Om» 

Uzz  -  - - t'mi- 


(3.70) 


The  horizontal  (x-  or  y-  directed)  dipole  source  may  give  rise  to  all  components  of  A  ,  i.e., 
A m>  =  xA™  +  yA ™  +  zA™‘.  The  transverse  electric  field  components  can  be  expressed  as 


d 


E™  =  (<  -  K)A™  -  kxkyA™  —  jkx—A 


(3.71) 
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Bmi 

•  =_m; 


-k,k„A?  +  (k2m  -  kl)A “  -  jkyU"- 


dz 


On  the  other  hand,  from  (3.20)  we  find 
E?'=jkxV^+jkyV*t 

K'=jkyV^-jkxV*t. 


(3-72) 

(3.73) 

(3.74) 


Eliminating  E™'  and  E™1  from  (3.71)-(3.74),  we  have  two  independent  equations  and  three 
unknown  components  of  A  .  Therefore,  the  solution  of  these  equations  is  not  unique,  and 
one  of  the  components  of  the  vector  potential  can  be  set  to  zero. 

3.3.1  Traditional  Form  of  the  Green’s  Function 

In  the  traditional  approach  [16],  it  is  postulated  that  an  x-directed  current  dipole  generates 
the  x  and  2  components  of  the  vector  potential.  Hence,  eliminating  E™  and  E™  from  (3.71)— 
(3.74),  and  letting  A™  =  0,  we  obtain 

kB 2 


■Vl 


Uky 


Ur  =  l 


k2  Ve _ -k2  Vh 

am  rmt  zm  mt 


dz  '  JUJ 

Substituting  (3.59)  and  (3.60)  into  (3.75)  and  (3.76),  we  have 


G~mi  _  1  riVh 

xx  ■  ^mt 


and 


hx 

fit2 

GVe  -  k 

dz  ZI 

ink  2 

[m 

mt  . 

which  we 

obtain 

nmi  —  - 

K 

\kl 

•— GVe 

^ zx  — 

U)kp 

k2 

L  zm 

dz°m' 

— GVh 

dz”1' 


(3.75) 

(3.76) 

(3.77) 

(3.78) 

(3.79) 


Similarly,  we  assume  that  the  y-directed  current  dipole  generates  the  y  and  z  components  of 
the  vector  potential.  Changing  x  to  y  and  kx  to  kv  in  (3.77)  and  (3.79),  we  obtain 


Gmt  jrivn  * 

yy  ^xx 


__  - 

'~7mi 


(3.80) 
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= _ K_  d  ^  _  d_^Vh 

zy  Ujkp2[k]mdz  m<  dz  m' 


(3.81) 


Thus,  the  elements  of  the  dyadic  (2.7)  have  been  obtained  in  the  Fourier  transform 
domain.  By  using  the  relationships  given  in  (3.8)-(3.13),  one  can  express  the  space  domain 
counterparts  of  (3.70)  and  (3.77)— (3.81 )  as 


/-rrru  /~imt 
^Jyy 


'  -  ?ks°  ^ 


/wit  _  1  _ c  j;  (  k m  ^  fiVc  _  ^  r^h 

°zx  ~  2tt  juj  4 51  [  kp 2  {  klm  dz  °mi  dz  G,n' 

G?v'  =  tan  C  G™J 


(3.82) 

(3.83) 

(3.84) 


G"'  =  2 l5"’]  • 

3.3.2  Alternative  Form  of  the  Green’s  Function 


(3.85) 


Taking  an  alternative  approach  [15],  we  suppose  that  a  horizontal  (x-  or  y-  directed)  dipole 
source  gives  rise  to  x  and  y  components  of  A™ .  Eliminating  E ™  and  E™  from  (3.7 1 )— (3.74), 
and  letting  A™  =  0,  we  obtain 


Af  =  -^~($-ik,VL-iKvL 

\Kzm 


(3.86) 

(3.87) 


If  the  voltages  V£i  and  V^,  in  (3.73)  and  (3.74)  are  those  generated  by  an  x-directed  dipole, 
then  A™'  =  G™  and  A™'  =  G™.  Hence,  substituting  from  (3.59)  and  (3.60)  into  (3.86) 
and  (3.87),  we  obtain 


Gm« 
yx 


L{ 

klkl 

ju  v 

ky 

J*L 

jukp 

2  Ul 

nVt  _i_  y  r*V s 

^mt  '  i  2  ^mt 

K  0 


^ mi 


(3.88) 

(3.89) 


Similarly,  substituting  from  (3.63)  and  (3.64)  into  (3.86)  and  (3.87),  we  obtain 


27 


(3.90) 


/orm  _  pmt 
^  xy  —  ^yx 


and 


1  /  P  P 

f~rm\  _  1  /  m  V  fiVc  i  "X  (~iVh  ) 

yy  J“\k2zmh2  mt  k/ Um'  y 


kl  ^ vh 


(3.91) 


That  G™  =  G™  was  already  anticipated  in  (2.8). 

Thus,  we  have  obtained  the  elements  of  (2.8)  in  the  Fourier  transform  domain.  Upon 
using  (3.79)— (3.13),  they  can  be  expressed  in  the  space  domain  as 


/imi  _  __2 _ 

xx  “  Airjw 


fr mi  ^ 

yy  ~  4  irju 


So 

So 


1.2 

Km  riVc  |  /iVh 
1,2  ^ mi  '  i 

L  Kzm 


L2 

Km  CV,  |  pVh 
j^2  '  “mi 


—  COS  ( S2 

. 

feM 

+  COS  (S 2 

v( 

kl  ~ 


Mm 


'  kl 

m 

KL 


vc 

mi 


We 

Tmi 


7  mi 


rvh 
'-7  mi 


l 


G”»'=G-=-4^8in2^ 
with  G™‘  given  by  (3.85). 


^  (  km  fiV'  _  p iVh 

k2  Ul  m> 


(3.92) 

(3.93) 

(3.94) 


3.4  Electric  Field  in  the  Layered  Medium  due  to  a 
Plane  Wave  Incident  in  Region  1 


Consider  an  incident  plane  wave  in  region  1 ,  given  as 
Ei(r)=(Eedi  +  E„<p1)e’k'-r 
where  the  propagation  vector  k\  is 

kx  =  k\  (sin$i  cos^  x  +  sin#i  sin^>j  y  +  cos#i  z) 


(3.95) 


(3.96) 


and  (<?i,<£j)  define  the  angle  of  arrival  of  the  plane  wave  in  terms  of  the  usual  spherical 
coordinate  convention.  We  can  utilize  the  analysis  in  Sections  3. 1-3. 3  to  obtain  the  total 
“incident”  electric  field  in  the  mth  layer  due  to  the  plane  wave  given  in  (3.95).  We  easily 
find  that  the  incident  field  in  the  top  layer  (m  =  1)  is  given  as 


Y*e(r)  =  {*  -£„sin^,  (l  + 

+  Eg  COS<^>i  COS  $1  (l  +  CO.*, (*-*,)' 
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+y  cos^  (l  + 

+£flsin*>,  cos0,  (l  +  r‘'e-j2fclCosfll(z“Jl)) 

—zEg  sin  fl  —  4p"ee_j2*:i  cosfil^_-Jl)Nj  1  ejil^in9l^C08v:’ir+sin'/>iy^",’c09^lJ^ 


and  the  field  in  other  layers  (m  /  1)  as 


EZ(r)  =  {*  [— ^ sin ipi  (l  +  T{)  T^z) 


+E,cos^-zt  l-  rn 

7]\  \  / 


+y  £„ cos v?i  (l  +  r*)  2^\(z)  +  ^sin^-j-Zt  (l  -  t\)  fv^(z) 


-zEe^  sin0m  (^1-  r?J  T^(2)|e 


jfcm(cos^i  9in9mx+sin^i  sin 9my+cos 9mz) 


(3.98) 


where  Z\ ,  F  1?  and  Tmi  are  given  in  Section  3.2,  with  kzm  =  kmcos9m,  where  9m  can  be 
obtained  from  the  Snell’s  law,  k\  sin$i  =  fcmsin#m.  In  (3.98).  rjm  =  \jjLmle-m  is  the  intrinsic 
impedance  of  the  dielectric  of  region  m. 
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Chapter  4 

Formulation  of  the  Mixed-Potential 
Electric  Field  Integral  Equation 


As  mentioned  in  Chapter  2,  the  decomposition  embodied  in  (2.13)  is  non-unique,  which 
means  that  an  infinite  number  of  different  MPIEs  is  possible.  We  observe,  however,  by 
referring  to  (2.16),  that  the  presence  of  the  “correction”  function  Pm'  has  the  undesirable 
effect  of  introducing  new  elements  in  the  dyadic  kernel  K™\  in  addition  to  those  already 
present  in  G™‘-  Ideally,  of  course,  Pmt  should  be  zero,  which  -  unfortunately  -  is  only 
possible  in  a  few  special  cases  [15],  Hence,  the  best  one  can  do  is  to  develop  MPIEs  for  which 
one  or  two  of  the  components  of  Pm'  are  zero.  It  is  shown  below  that  for  the  vector  potential 
Green’s  functions  (2.7)  and  (2.8)  the  x  and  y  components  of  Pmt  are  not  independent,  thus 
leaving  us,  for  each  G™' ,  with  only  two  degrees  of  freedom:  either  P™'  =  P™1  =  0  and 
Pzm‘  ^  0,  or  P™'  0,  P™'  ^  0,  and  P™  =  0.  Of  the  four  possible  “minimal”  formulations, 

only  three  lead  to  distinct  MPIEs.  These  three  are  discussed  in  detail  below,  where  they  are 
referred  to  as  Formulations  A,  B,  and  C.  Most  of  our  derivations  are  done  in  the  Fourier 
transform  domain  defined  in  (3.1)  and  (3.2),  which  greatly  simplifies  the  algebra. 

4.1  Formulation  A 

In  this  formulation,  we  employ  the  alternative  form  of  G™  given  in  (2.8).  In  the  Fourier 
transform  domain,  the  x ,  y,  and  z  components  of  (2.13)  become 


w  =  hk,K™  f 


(-jk.ee;  =  —jk,i<r  +  juP; 
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(4.3) 


JU,  3  xm<-  _  4  9  tsttu  -  •  nrai 

W dlG-  -]ZaT'ht  +,u>p-  ■ 


Using  (3.88)  through  (3.91)  in  (4.1)  and  (4.2),  one  finds  that  P™'  and  Pymt  are  related  by 

(4.4) 


~  k  ~ 

pmi  _  V  pmi 

y  L  * 

rK  7* 


In  Formulation  A,  we  choose  P™‘  =  P™  =  0,  in  which  case  A'”'  can  be  interpreted  as 
the  scalar  potential  of  a  point  charge  associated  with  a  horizontal  dipole  [I5j.  Solving  (4.1) 
or  (4.2)  for  A'™1,  one  obtains 

GVJ, 


A"77U 

0 


k* 


which  can  be  substituted  into  (4.3)  to  yield 

pmi  _  /*■<-!  ~  Fm  (-m  ry; 
etKzm 

where  we  have  introduced  for  later  convenience 


(4.5) 


(4.6) 


= 

■£m  i 


zl-iGu  i  -  1  >  m  >  1 


(4.7) 


-  ZgtG^(zt,  z')  7^(2),  n  +  1  >  m  >  1  +  1 

in  which  the  superscript  “<?”  stands  for  “e”  or  “/i.”  Observe  that  =  0  when  m  =  i,  i.e.. 
when  the  source  and  observation  points  are  within  the  same  layer. 

Substituting  (3.70),  (3.88)— (3.91 )  and  (4.6)  into  the  Fourier  domain  counterpart  of  (2.16), 
and  using  the  relations  given  in  (3.8)— (3.13),  one  obtains  the  dyadic  kernel 

K™(r  \r')  =  x  xK™  +  y  yK™  +  »  zK™ 

+(*  V  +  ifx) K™  +  *  zK™  +  y  zK ™  (4.8) 

with  the  elements  given  by 


A" = gs  - {*is»  (4:<5”‘) + s»  os) 


+  cos  2(  S2 


—  ( rV»  _  P'V'  " 
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(4.10) 
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A«=G‘”'  =  4^sin2<Sj 
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l  kp 
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at;  =  S- p™  =  -^fs1  (Pr) 


dx 
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Kmt  —  __pmi  — 
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2  7T 

sin  ( 
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5i  (AT) 
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-  G«  +  =  27nWr 


■5o  (GV)  • 


(4.11) 

(4.12) 

(4.13) 

(4.14) 


Finally,  the  Fourier  inversion  of  K™  in  (4.5)  yields  the  scalar  potential  kernel 


(4.15) 


We  observe  from  (4.9)  through  (4.15)  that  in  this  formulation,  when  m  ^  i,  the  effect  of 
V  Pm'  in  (2.16)  is  to  introduce  two  new  entries,  K™  and  AT\  and  to  modify  G™ .  However, 
when  the  object  is  confined  to  a  single  layer  (m  =  1),  we  simply  have  K =  G™1,  so  no 
modification  of  the  Green’s  function  is  required. 

A  useful  property  of  Formulation  A  is  the  cancellation  of  the  contour  integrals  in  (2.15), 
which  is  the  result  of  the  continuity  with  respect  to  the  z'  coordinate  of  the  scalar  potential 
kernel  at  the  ith  interface:  K™l{z'  =  Zi  -f  0)  =  K™',+1(z'  =  z,  —  0).  We  note,  however, 
that  a  continuity  with  respect  to  the  z  coordinate  does  not  hold,  he.,  K™(z  =  zm  +  0)  7^ 

A7+1,i(z  =  zm  -  0). 


4.2  Formulation  B 


In  this  formulation,  as  in  Formulation  A,  we  employ  the  alternative  form  (2.8)  of  G™'- 
However,  rather  than  choosing  P™  =  P™’  =  0,  we  select  P™1  =  0  in  (4.1)  through  (4.3).  In 
Inis  case  K ™  can  be  interpreted  as  the  scalar  potential  of  a  point  charge  associated  with  a 
vertical  dipole  [15].  From  (4.3)  and  (3.70),  we  find 


=  -j« 


.  G 


K 


*  JW  kl 

which  can  be  substituted  into  (4.1)  to  yield 


nm>  _  1.  MmeT7i  Piei  ^Ve 

rx  -  ljkx  umi 

KziKzm 


(4.16) 


(4.17) 
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which,  in  view  of  (4.4),  also  specifies  P™'. 

Referring  to  (2.16)  and  proceeding  as  in  Formulation  A,  we  obtain  the  dyadic  kernel 


KT(r  \r')  =  x  xK™  +  y  yK™  +  x  zK 
+  {xy  +  yx)K™  +  zxK™  +  zyKi 

with  the  elements  given  by 
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where  P7"'  is  nonzero  only  for  m  ^  i  and  is  given  by 

nmi  _  P«£«  ~  /*m£ro  r/, 

n  e  k2  mi 

cm,v2  f 


(4.20) 

(4.21) 

(4.22) 

(4.23) 

(4.24) 


(4.25) 


in  which  7^  can  be  obtained  from  7^  given  in  (4.7)  by  replacing  in  the  latter  the  charac¬ 
teristic  impedances  Z\  by  their  reciprocals. 

Finally,  the  Fourier  inversion  of  K™  in  (4.16)  gives  the  scalar  potential  kernel 


jy-mt  3^  o  (  ^rni 

K*  --2^°Ur 


(4.26) 
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We  note  that  in  Formulation  B  two  new  entries  (K™  and  iv™‘),  not  present  in 


imi 

Va  ' 


are  introduced  to  the  dyadic  kernel  K™' .  Also,  extra  terms  are  added  to  G™,  G™\  and 


yy ' 

G™'.  In  the  case  where  the  object  is  confined  to  a  single  layer  (m  =  i),  K™  =  G™‘ ,  and 
Formulation  B  becomes  identical  to  Formulation  A. 

The  continuity  properties  of  the  scalar  potential  kernel  in  Formulation  B  are  complemen¬ 
tary  to  those  in  Formulation  A;  that  is,  in  the  present  case  K™'(z'  =  z,  +  0)  ^  A'™,1+1(z'  = 
z,  —  0)  and  K™(z  =  zm  +  0)  =  K™+1,,(z  =  zm  —  0).  As  a  result,  the  contour  integrals 
in  (2.15)  do  exist  when  the  object  penetrates  one  or  more  of  the  interfaces. 

4.3  Formulation  C 

The  point  of  departure  in  this  formulation  is  the  traditional  form  of  G™1  given  in  (2.7). 
Using  it  in  (2.13),  we  obtain  in  the  Fourier  transform  domain 


• 

jw  / 

ki  l 

ju  ( 

kl  V 

• 

vy'“rxr 
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+ jup: 


(4.27) 

(4.28) 

(4.29) 
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From  (4.27)  and  (4.29),  and  referring  to  (3.79)  through  (3.81),  we  find  that  (4.4)  still  holds. 
Since  G™'  in  the  present  case  is  the  same  as  that  in  the  “alternative”  Green’s  function  (2.8), 
specifying  =  0  results  in  an  MPIE  identical  to  that  in  Formulation  B.  We  therefore  set 
=  0  in  (4.27)  and  (4.28),  which  yields  the  scalar  potential  kernel 

(<5£  -  61 V)  .  (4.30) 

From  the  above  and  (4.29),  there  results 

+  (4-31) 

Substituting  (3.70),  (3.79)  through  (3.81)  and  (4.31)  into  the  Fourier  domain  counterpart 
of  (2.16),  and  using  the  relations  given  in  (3.8)— (3.13),  we  obtain  the  dyadic  kernel 

K7(r  |r')  =  {xx  +  yy)  K™  +  x  zK™  +  y  zK™ 
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(4.32) 


with  the  elements  given  by 
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(4.34) 

(4.35) 
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with  1%.  given  in  (4.7).  For  notational  simplicity,  we  have  suppressed  in  (4.39)  and  (4.40) 
the  superscripts  “e”  or  “h”.  We  note  that  one  can  obtain  WP  from  by  replacing  in  the 
latter  the  characteristic  impedances  Z,  by  their  reciprocals. 

Finally,  the  Fourier  inversion  of  (4.30)  gives 


l/mi  7^  q 
K *  =  2^° 


A  (%  -  ea) 


(4.41) 


As  in  Formulation  A,  this  K™  may  be  interpreted  as  the  scalar  potential  of  a  point  charge 
associated  with  a  horizontal  dipole  [15].  However,  these  two  potentials  are  not  identical, 
since  each  corresponds  to  a  different  form  of  the  vector  potential  Green’s  function. 

We  observe  from  (4.33)  through  (4.38)  that  Formulation  C  introduces  two  new  entries 
(A™‘  and  K ™),  not  present  in  Q.™,  to  the  dyadic  kernel.  Also,  extra  terms  are  added  to 
G™1.  In  contrast  to  Formulations  A  and  B,  these  modifications  occur  even  if  the  object  is 
confined  to  a  single  layer.  As  in  Formulation  A,  the  scalar  potential  kernel  in  the  present 
case  has  the  continuity  property  that  K™‘(z'  =  zx  -f  0)  =  K™'l+l(z'  =  2,  —  0),  which  causes 
the  contour  integrals  in  (2.15)  to  cancel.  Formulation  C  also  shares  with  Formulation  B  the 
useful  property  that  K™'(z  =  zm  +  0)  =  K™+1,'(z  =  zm  —  0). 
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4.4  Discussion 


The  properties  of  three  MPIEs  are  summarized  for  easy  reference  in  Table  4.1.  Each  of 
the  three  formulations,  called  A,  B  or  C,  can  be  derived  from  either  the  alternative  or  the 
traditional  form  of  the  dyadic  Green’s  function  for  the  vector  potential.  These  two  dyadics 
are  shown  in  matrix  form  in  the  first  column  of  Table  4.1.  We  note  that  Formulation  A  is 
derived  from  the  alternative  forms  of  Formulation  C  from  the  traditional  form,  and 
Formulation  B  from  either  of  the  two  (this  is  due  to  the  fact  that  both  forms  of  share 
the  same  G™').  The  distinguishing  feature  of  each  of  the  three  formulations  is  the  choice  of 
the  scalar  potential  kernel  K™',  which  also  specifies  the  “correction”  vector  Pmx  according 
to  (2.13).  Although,  as  the  second  column  of  Table  4.1  indicates,  the  scalar  potential  kernels 
in  Formulations  A  and  C  are  both  associated  with  a  horizontal  dipole,  they  are  different, 
because  they  correspond  to  different  vector  potential  Green’s  functions.  In  Formulation  B, 
K™1  is  that  associated  with  a  vertical  dipole.  Actually,  by  properly  choosing  K?',  one  can 
also  derive  Formulation  A  from  the  traditional  form  of  G™',  and  Formulation  C  from  the 
alternative  form. 

The  forms  of  the  dyadic  vector  potential  kernel  are  shown  for  each  of  the  three  formula¬ 
tions  in  the  third  column  of  Table  4.1.  Comparing  this  column  with  column  one,  we  observe 
that  in  all  three  formulations  two  new  entries,  in  addition  to  those  already  present  in 
appear  in  K ™‘,  which  is  of  course  undesirable.  We  note,  however  (as  indicated  in  the  fourth 
column  of  Table  4.1),  that  in  Formulations  A  and  B  the  number  of  the  entries  in  K™'  is 
not  increased  over  that  in  when  the  object  is  confined  to  a  single  layer  (m  =  i ).  We 
should  also  point  out  that  the  correspondence  between  the  number  of  entries  in  K™'  and 
the  number  of  distinct  Sommerfeld  integrals  that  need  be  evaluated  is  not  one-to-one.  In 
fact,  we  can  show  by  referring  to  Sections  4. 1-4. 3  that  in  the  general  case  only  four  distinct 
Sommerfeld  integrals  are  required  in  all  three  formulations.  When  m  =  i,  the  number  of 
distinct  integrals  in  Formulations  A  and  B  reduces  to  three. 

As  mentioned  at  the  end  of  Chapter  2,  one  may  leave  the  Pm‘  term  as  a  part  of  the 
scalar  potential  (cf.  (2.15)),  thus  avoiding  the  modification  of  the  vector  potential  kernel  as 
in  (2.16).  However,  since  this  would  constitute  a  departure  from  the  standard  form  of  the 
MPIE  [4,  5],  we  prefer  to  proceed  according  to  (2.16).  Although  in  our  formulations  some 
of  the  terms  introduced  to  K™  by  V  Pmi  in  (2.16)  become  singular  when  r  and  r’  coincide 
on  an  interface,  they  are  no  more  singular  than  the  terms  already  present  in  and  can 
be  handled  in  a  similar  way  as  the  latter. 
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The  continuity  properties  of  the  scalar  potential  kernels  across  the  interfaces  are  sum¬ 
marized  in  the  fifth  column  of  Table  4.1.  In  Formulations  A  and  C,  K™'  is  continuous  with 
respect  to  z1,  which  results  in  the  cancellation  of  the  undesirable  contour  integrals  in  (2.15) 
(c/.the  sixth  column).  In  Formulations  B  and  C,  K™'  is  continuous  with  respect  to  z,  which 
results  in  considerable  simplifications  in  the  numerical  procedure  when  the  object  penetrates 
one  or  more  interfaces.  This  last  point  can  be  fully  appreciated  only  after  we  have  discussed 
the  details  of  the  MOM  procedure  in  the  next  chapter. 

We  conclude  from  the  above  summary  that  when  the  object  is  confined  to  a  single  layer, 
Formulations  A  and  B  become  identical  and  are  preferable  to  Formulation  C,  because  they 
have  fewer  terms  in  K ™.  In  the  general  case,  Formulation  C  enjoys  a  clear  advantage 
over  Formulations  A  and  B,  because  it  does  not  have  contour  integrals  and  because  its 
scalar  potential  kernel  is  continuous  at  the  interfaces  with  respect  to  2,  which  results  in 
the  simplifications  mentioned  above.  If  we  had  to  choose  between  Formulations  A  and  B, 
we  would  prefer  the  latter,  because  the  advantages  of  having  the  scalar  potential  kernel 
continuous  with  respect  to  2  more  than  compensate  for  the  complications  caused  by  the 
presence  of  the  contour  integrals.  This  point  is  further  elaborated  upon  in  the  next  chapter. 

The  previous  works  related  to  the  mixed-potential  EFIE  and  reviewed  in  Chapter  1  can 
be  classified  as  shown  in  the  last  column  of  Table  4.1. 
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Chapter  5 

Numerical  Method 


In  the  previous  chapter,  three  alternative  MPIE  formulations  -  referred  to  as  Formulations 
A,  B  and  C  -  have  been  derived  for  PEC  antennas  or  scatterers  residing  in  plane-stratifieu 
dielectric  media  with  an  arbitrary  number  of  layers.  Formulation  C  has  been  found  to  be 
particularly  well  suited  for  the  application  of  the  MOM  procedures  originated  by  Wilton 
and  his  co-workers  [4,  5,  7].  This  formulation  is  adopted  here  in  conjunction  with  those 
procedures  to  analyze  the  radiation  and  scattering  by  a  PEC  object  of  arbitrary  shape 
residing  in  layered  media.  For  simplicity,  but  without  loss  of  generality,  we  describe  this 
method  in  the  context  of  the  two-media  problem  illustrated  in  Fig.  5.1.  Also,  for  simplicity, 
the  interface  between  the  two  media  is  taken  to  be  the  ary-plane  of  a  Cartesian  coordinate 
system.  The  upper  (z  >  0)  half-space,  which  is  characterized  by  t\  and  /j1?  will  be  referred 
to  as  region  1,  and  the  lower  ( z  <  0)  half  space  in  Fig.  5.1a,  or  the  slab  (—d  <  z  <  0)  in 
Fig.  5.16,  both  characterized  by  e2  and  /z2,  as  region  2.  The  parts  of  the  surface  S  of  the 
object  which  are  in  regions  1  and  2  are  denoted  as  S\  and  52,  respectively. 

As  was  already  mentioned,  based  on  the  analysis  of  Section  4.4,  Formulation  C  has  been 
selected  for  the  problems  of  Fig.  5.1.  Hence,  upon  specializing  the  MPIE  of  (2.15)  to  the 
present  case,  one  obtains 
2 

n0x^[jwr(r)+Vr(r)]  =  hm  x  E™  (r) ,  r  on  Sm,  m  =  l,2  (5.1) 

t=i 

where,  for  notational  simplicity,  we  have  introduced 

Ami(T)  =  jjQi{r\r')-J(r')dS’  (5.2) 

s. 

r‘(')  =  j  K?(,T\r')q(r')dS'  (5.3) 

s, 
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/  /  /  / 


Figure  5.1:  PEC  object  of  arbitrary  shape  embedded  in  (a)  contiguous  half-spaces,  and  (b) 
in  a  grounded  slab. 


in  which  the  kernel  elements  are  given  in  Appendix  A.  The  integrals  (5.2)  and  (5.3)  represent 

#  the  modified  potentials,  and  should  not  be  confused  with  the  original  potentials  given  in  (2.4) 
and  (2.5). 

We  note  that  the  contour  integrals  encountered  in  (2.15)  are  absent  in  (5.1),  which  is 
characteristic  of  hormulation  C.  The  quantity  of  interest  for  the  scattering  and  radiation 

#  problems  is  the  surface  current  distribution  on  S,  from  which  other  characteristics,  such  as 
the  far  field  pattern  or  radar  cross-section,  can  readily  be  computed.  For  the  transmission 
line  problem,  there  is  no  excitation  and  (5.1)  becomes  a  non-standard  eigenvalue  problem, 
with  the  propagation  constants  as  eigenvalues  and  the  modal  current  distributions  as  eigen- 

•  functions. 

In  solving  the  integral  equations,  the  MOM  schemes  developed  in  [3,  4,  5,  7]  are  employed. 
These  procedures  are  modified  to  account  for  the  dyadic  character  of  the  vector  potential 
kernel,  to  ensure  the  current  continuity  across  the  interface,  and  to  allow  for  the  charge 

•  discontinuity  there  [19].  In  the  MOM,  basis  functions. are  chosen  to  represent  the  unknown 
currents,  testing  functions  are  chosen  to  enforce  the  integral  equation,  and  these  are  used  to 
derive  an  impedance  matrix  approximant  to  the  integral  equation. 


5.1  Surface  of  Arbitrary  Shape 

Following  the  procedure  described  by  Rao  et  al.  [5],  the  surface  S  of  the  ob,  ~t  is  modeled  in 
terms  of  triangular  patches  in  a  manner  suggested  in  Fig.  5.2.  The  surface  current  density 
on  S  is  approximated  as 


N 


J(r)  =  £  /,A.( r) 


(5.4) 


n=l 


where  N  is  the  number  of  interior  (nonboundary)  edges  and  An  is  the  vector  basis  function 
defined  on  the  adjacent  triangles  associated  with  the  nth  edge,  and  is  given  as 


Ai(r)  = 


—  r  in  S+ 
*+’  " 

7~ *  T  in  S~ 

hn 

(  0,  otherwise 


(5.5) 


in  which  h±  is  the  height  of  triangle  S*  relative  to  the  nth  edge  of  length  and  ±p„  is  the 
vector  from  the  free  vertex  of  S„  to  an  arbitrary  point  on  patch  ,  as  shown  in  Fig.  5.3.  The 
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Figure  5.3:  Local  coordinates  associated  with  an  e 


expansion  for  the  surface  charge  densitj  q ,  which  is  obtained  by  using  (5.4)  in  the  equation 
of  continuity,  calls  for  the  surface  divergence  of  An,  which  is  found  as 


V,-An(r] 


r  in  S+ 

r  in  S~ 
otherwise  . 


(5.6) 


When  the  expansions  for  J  and  q  are  used  in  (5.1)  and  the  resulting  equations  tested  by 
integration  along  the  path  from  the  centroid  of  S+  to  the  middle  of  the  edge  lp  and  then  to 
G  the  centroid  of  S~,  an  N  x  N  system  of  linear  equations  is  obtained,  which  may  be  written 

in  matrix  form  as 


[ZJ[/J  =  K1 

(5.7) 

• 

where  the  elements  of  Zpn  and  Vp 

are  given  as 

ZPn  -  y  (Apn  +  Ap„)  -  —  ($pn  -  $+„) 

(5.8) 

• 

K  =  J  (Er-  +  E+) 

(5.9) 

where  p,  n  =  1 , 2, . . . ,  N,  and 

• 

A%  =  JpfK 7  (r?|r 

S, 

')  •  An  (r’)  dS1 

(5.10) 

9% = j  K' « i  '•') 

S, 

•  An  (r')  dS’ 

(5.11) 

# 

Et  =  P?  •  EZ  (rf)  ■ 

(5.12) 

In  the  above,  pc^  is  the  vector  between  the  free  vertex  and  the  centroid  of  Sp,  with  pp~ 
directed  toward  and  pcp+  directed  away  from  the  vertex,  and  rc±  is  the  vector  from  the 
coordinate  origin  to  the  centroid  of  Sp  . 

It  is  worthwhile  to  elaborate  on  one  aspect  of  the  above  procedure,  which  is  peculiar  to 
the  layered  medium  case.  As  suggested  in  Fig.  5.2,  the  triangle-patch  model  of  5  must  take 
into  account  the  position  of  the  interface,  since  no  triangular  patch  is  allowed  to  cross  it. 
This  arrangement  ensures  the  continuity  of  the  current  on  S  across  the  interface  and  allows 
for  the  discontinuity  of  the  charge  there  (cf.  (5.5),  (5.6)).  A  typical  triangle  pair  associated 
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with  edge  p  on  the  interface  is  illustrated  in  Fig.  5.4.  In  testing  (5.1)  with  Ap ,  one  should 
bear  in  mind  that  S~  and  S£  are  in  different  layers.  For  example,  testing  Vd>mi  in  (5.1) 
with  Ap  yields 

J  V<t>2'  ■  ApdS  +  j  V<j>u  ■  ApdS  =  -  j  4>2,V  ■  Ap  dS  -  J  <t>u  V  •  Ap  dS 

Sp  Sp  Sp  Sp 

+  /  (d>%+  •  Up  +  <f>uA~  ■  ii~ )  dl  (5.13) 

where  use  has  been  made  of  the  Gauss’  theorem  [38,  p.  503]  and  where  Ap  and  A~  signify, 
respectively,  Ap  evaluated  at  the  top  and  bottom  sides  of  the  interface.  We  observe  the 
undesirable  line  integral  in  (5.13),  which  does  not  occur  for  objects  in  homogeneous  space. 
Even  in  the  present  case,  however,  this  integral  disappears  for  Formulations  B  and  C  (see 
Chapter  4)  in  view  of  the  continuity  of  the  scalar  potential  kernels  as  e  crosses  an  interface 
and  the  fact  that  ■  A+  =  —  u~  •  A~  on  lp,  which  follows  from  (5.5).  In  Formulation  A, 
the  line  integrals  contribute  a  contour  integral  at  the  intersection  of  5  with  the  interface, 
quite  similar  to  that  in  (2.14),  which  is  present  in  Formulation  B.  However,  due  to  the 
approximations  made  in  the  testing  procedure  [5],  the  latter  is  actually  easier  to  handle  than 
the  former-hence  our  statement  in  Chapter  4  regarding  the  superiority  of  Formulation  B 
over  A.  Formulation  C  is,  of  course,  free  of  either  of  these  contour  integrals. 

It  is  worth  noting  that,  apart  from  the  presence  of  the  Sommerfeld  integrals,  the  only 
difference  between  (5.8)— (5.12)  and  their  free  space  counterparts  [5]  is  the  dyadic  character 
of  the  vector  potential  kernel  in  (5.10).  The  surface  integrals  in  (5.10)— (5.1 1 )  can  be  effi¬ 
ciently  evaluated  in  normalized  area  coordinates  and  reused  for  different  matrix  elements,  as 
described  in  [5]. 

As  in  the  free  space  problem,  the  integrands  of  Am‘(r)  and  $m'(r)  are  singular  at  r  =  r', 
i.e.,  when  the  observation  point  coincides  with  the  source  point.  As  mentioned  in  Ap¬ 
pendix  A,  this  singularity  arises  from  the  closed  form  part  of  the  Green’s  function.  Thus, 
when  the  observation  point  is  close  to  the  source  point,  extraction  of  the  inverse-|  r  —  r'  | 
behavior  is  necessary  before  numerical  integration  is  performed.  The  technique  developed  in 
[41]  may  be  used  for  the  treatment  of  such  singularities. 

5.2  Thin  Wire  of  Arbitrary  Shape 

We  next  consider  a  thin-wire  structure  of  arbitrary  shape  with  radius  a,  which  may  penetrate 
the  interface  between  two  dissimilar  dielectric  layers,  as  shown  in  Fig.  5.5.  The  parts  of  the 
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2  a 


JV  +  1 


Figure  5.5:  Straight-segment  model  of  a  thin-wire  of  arbitrary  shape  which  penetrates  the 
interface  between  dissimilar  media. 
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wire  which  are  in  regions  1  and  2  are  denoted  as  L\  and  £2,  respectively.  The  wire  is 
modeled  in  terms  of  N  +  1  straight  tubular  segments  and  the  vectors  r„  n  =  0,  1, 

N  +  1,  defined  with  respect  to  the  global  coordinate  origin  O,  locate  the  end  points  of  the 
linear  segments,  as  illustrated  in  Fig.  5.5.  The  same  points  are  designated  by  /„  with  respect 
to  the  local  arc-length  coordinate  /  along  the  axis  of  the  wire.  By  making  use  of  the  thin 
wire  approximation  [3],  we  may  express  the  potentials  (5.2)  and  (5.3)  as 


Amt(i)  —  j  k™(i  1  /')•/(  i')  dr 

L, 

4>mi(/)  =  J  K?'(i\i')Q(r)dr 


(5.14) 


(5.15) 


where  we  have  introduced  the  total  current  1(1)  =  2iraJ(l)  and  the  linear  charge  density 
Q(l)  =  2x aq(l).  Since  on  the  wire  r  is  determined  by  /,  we  have,  for  simplicity,  replaced  the 
former  by  the  latter  in  (5.14)  and  (5.15).  Similar  remarks  apply  to  r'  and  Following  the 
MOM  developed  in  [7],  the  current  density  along  the  wire  is  approximated  as 


N 


1(1)  =  £  InAn(l) 

n=  1 

where  An(l)  is  the  vector  expansion  function  given  by 


(5.16) 


p t 

/  in  5+ 

• 

Ml)  =  < 

Pn 

l  in  5~ 

(5.17) 

.  0, 

otherwise  . 

in  which  h±  is  the  length  of  the  segment  5±  relative  to  the  nth  node  and  ±p*  is  the  vector 
from  the  free  node  of  to  an  arbitrary  point  on  the  segment,  as  shown  in  Fig.  5.6.  The 
divergence  of  An(l),  which  is  proportional  to  the  linear  charge  density  assc~;ated  with  this 
basis  function,  is 


V  •  An(l)  = 


K' 
1 

0, 


l  in  5+ 


—  ,  /  in  Si 


(5.18) 


otherwise 
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nth  node 


Figure  5.6:  Local  coordinates  associated  with  the  nth  node  of 
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When  the  expansions  for  I  and  Q  are  used  in  (5.1)  and  the  resulting  equations  tested  with 
A„(l),  an  N  x  N  system  of  linear  equations  is  obtained,  which  may  be  written  in  matrix 
form  as 


=  ra 


where 


[K  +  a-„)  -  i-  (*-  -  *+) 


vP  =  Pi+-E':v;+)+pr-ET(n-) 


in  which 


A%=  J  pf  '  10-  ^40  dl’ 

L, 

*£  =  /  10  v'-An(i')dr 


(5.19) 

(5.20) 

(5.21) 

(5.22) 

(5.23) 


and 


_  Ip 

p  2 


(5.24) 


In  (5.22)  and  (5.23),  we  employ  the  exact  kernel  on  the  source  segment  and  the  reduced 
kernel  otherwise. 

It  is  observed  that  the  integrand  of  (5.22)  consists  of  seven  terms  contributed  by  the 
non-zero  elements  of  the  dyadic  kernel  For  example,  the  contribution  of  K™  to  A ^  is 

/  p?  ■  xKZV?  1 0  z  ■  M 0  dl’ 

l n  ^n+1 

/  K™(i?\o^r-zdr+  J  K™(i;+\r)*t-zdr 


=  pcv  ■  * 

fcl.  i 

The  integrals  in  (5.25)  may  be  further  transformed  as 

n'+ 


(5.25) 


In 


J  K™(l?\l')^dl' 

l-l 

A+  A+ 

=  2  /  Kr;(';+i»'+o*'  +  ^  /  +£>'*' 


(5.26) 


-A+ 


-A+ 
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and 


*n+  1  f_ 

In  U 

-  A"  A- 


(5.27) 


where  A±  =  h^/2. 

An  examination  of  (5.20)  and  of  the  above  expressions  reveals  that  the  right-hand  inte¬ 
grals  in  (5.26)  can  be  used  in  the  matrix  elements  Zp_1)n_i,  ZPi„_ i,  Zp_i,n,  and  Zpn,  resulting 
in  considerable  savings  in  matrix  fill  time. 

As  suggested  in  Fig.  5.5,  no  straight  segment  is  allowed  to  cross  the  interface,  where  the 
charge  is  discontinuous.  Also,  by  using  the  same  arguments  as  given  in  Section  5.1,  one 
can  prove  that  no  extra  terms  need  to  be  added  to  (5.20)  to  account  for  the  interface,  if 
Formulation  C  is  used. 

5.3  Microstrip  Patch  Antenna  of  Arbitrary  Shape 

The  structure  considered  in  this  section  is  a  microstrip  patch  antenna  of  arbitrary  shape 
driven  by  a  coaxial  transmission  line,  as  illustrated  in  Fig.  5.7.  The  dielectric  slab  and  the 
PEC  ground  plane  are  assumed  to  be  of  infinite  extent.  The  shape  of  the  patch  (also  assumed 
to  be  PEC)  is  arbitrary  and  is  modeled  by  triangular  elements.  The  probe-to-patch  junction 
can  be  located  anywhere  on  the  patch,  including  edges  and  corners,  but  it  must  coincide 
with  a  node  of  the  triangle-element  mesh. 

Following  the  procedure  developed  by  Hwu  and  Wilton  [7],  the  surface  current  density  on 
the  patch  Sp  and  the  total  axial  current  on  the  wire  Sw  are  approximated  as,  respectively, 


J(r)  »  £  IPAp(r)  +  IJAJ(r),  r  on  SP 


(5.28) 


I(r)l  *  £  I*Aw(r)  +  IJAJ(r),  r  on  Sw 


(5.29) 


where  Ap(r)  and  Aw(r)  are  the  basis  functions  previously  defined  in  (5.5)  and  (5.17), 
respectively  (the  superscripts  “P”  and  “W”  have  been  introduced  to  distinguish  between 
functions  associated  with  the  patch  and  the  wire).  The  basis  function  associated  with  the 
wire-to-patch  junction  is  given  as  [7,  42] 
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microstrip  patch 


Figure  5.7:  Coax-fed  microstrip  patch  antenna,  (a)  Cross-sectional  view,  (b)  Top  view 
showing  a  triangle-element  model  of  an  arbitrarily  shaped  patch. 
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ylV)  =  { 


K, 


{hJn2 


i  - 


(pJr-h7)2 


i  J-  , 


<(r), 
l  0, 


Af(r),  r  in  Sf 

r  in  SJ+ 
otherwise 


(5.30) 


where  the  index  /  refers  to  the  /th  triangle,  Sf  ,  at  the  junction  (cf.  Fig.  5.8),  Af(r)  and 

A  J- 

h[  are  the  patch  basis  function  and  the  vector  height,  respectively,  associated  with  the 
edge  opposite  the  junction  vertex  in  Sf~,  A™ (r)  is  the  wire  basis  function  associated  with 
the  junction  node  Nw  +  1,  and  Kt  is  the  normalization  constant  given  as 


K, 


on 

Nj 

i=i 


on 

tlQ.1 


(5.31) 


where  cti  is  the  angle  between  the  two  edges  of  Sf~  common  to  the  junction  vertex,  /;  is  the 
length  of  the  edge  opposite  the  junction  in  Sf~,  Nj  is  the  number  of  triangular  elements 
attached  to  the  junction,  and  a 1  is  the  sum  of  the  Nj  junction  vertex  angles.  The  divergence 
of  AJ(r),  which  is  proportional  to  the  surface  charge  density  near  the  junction,  is  given  as 


V,  •  AJ{r) 


2  Kx 

1 

hJ+' 


r  in  Sf 
r  in  SJ+ 


(5.32) 


[  0,  otherwise. 


The  basis  function.  (5.30)  enforces  the  current  continuity  at  the  junction  and  it  correctly 
models  the  near-singular  behavior  of  the  patch  current  near  the  feed  probe.  Unlike  some 
other  “attachment  modes”  described  in  the  literature  [18,  43,  44,  45,  46],  it  is  applicable 
even  in  cases  where  the  microstrip  patch  antenna  is  driven  at  an  edge  or  a  vertex. 

The  testing  procedures  for  Sp  and  Sw  are  identical  to  those  described  in  Sections  5.1 
and  5.2,  respectively.  In  testing  the  equations  associated  with  the  junction,  we  first  integrate 
along  the  wire  axis  from  the  center  of  the  attached  wire  segment  to  the  junction,  and  then 
along  the  paths  from  the  junction  to  the  centroid  of  each  triangle  element  associated  with 
the  junction.  The  resulting  equations  are  subsequently  combined  into  a  single  equation  for 
the  junction,  by  weighting  each  with  the  associated  triangle  vertex  angle  and  summing  the 
results  for  each  junction  triangle  element.  Substituting  the  current  expansions  (5.28-5.29) 
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into  (5.1)  and  testing  the  resulting  equations  by  the  method  described  above,  we  obtain  an 
N  x  N  system  of  linear  algebraic  equations,  where  N  =  Np  +  Nw  +  1 ,  which  may  be  written 
in  matrix  form  as 


■  [zpp]  [zpw' 

M  ■ 

|^p]  gWW 

[zWJ  j 

.  [zJP\  [zJW' 

ZJJ 

Ip 

jW 

IJ 


\vw 


VJ 


(5.33) 


where  the  elements  of  the  submatrices  [Zpp  ,  ZWWJ,  [Vpj ,  and  |Vlv]  are  defined  in  pre¬ 
vious  sections,  and  the  elements  of  the  other  submatrices  are  given  as 


C  =  y  (Ap-  +  Ap+)  -  j ■-  (*p-  -  ip*)  ,  7  =  w  or  J  (5.34) 

Zp  =  (Ap*  +  Ap-)  -  y  (ip-  -  ip*)  ,  7  =  P  or  J  (5.35) 

J  ^ 


Zp  =  a,  (jwAfp  -  —ifr)  +  ju,Ap*  +  -ip*,  7  =  P,  W  or  J  (5.36) 

a  1= 1  V  J 

vJ  =  jy  pP  •  £"“(■•?-)  +  pcnw+\  •  £in‘('XU)  (5-37) 


in  which 


Ap*  =  /  P?  •  I  *•')  •  4!(0  1  =  W  or  J 

s 

Ap*  =  /  pf  •  £?((?  I  <■')  •  AIM  ds',  T  =  Pot  J 
s 

AJ~t-  _  _  4P'7_  4  J"r+  _  ^ W-r-j- 

'*/n  ^/n  ’  '‘Nw  +  l.ii 

*5*  =  /  A-f  (r  f  I  r')  V'  •  ^J(r')  7  =  (V  or  J 

S 

ip*  =  /  Kf(lf  |  r')  V'  •  AI(t')  7  =  P  or  J 


$ 


•/Tr— 

/n 


=  $f'y_ 

^(n  1 


=  $ 


W-r+ 
/Vw  +  l,n 


(5.38) 

(5.39) 

(5.40) 

(5.41) 

(5.42) 


(5.43) 


The  antenna  is  excited  by  a  coaxial  cable  from  below  the  ground  plane  (cf. Fig.  5.7). 
The  magnetic  current  frill  model  [47]  is  .  -~d  to  model  the  junction  between  the  coaxial 
transmission  line  and  the  dielectric  substrate.  It  is  assumed  that  the  field  distribution  in  the 
coax  aperture  is  that  of  the  TEM  mode  with  a  known  voltage  Vin.  The  magnetic  current 
radiating  in  the  grounded  slab' environment  provides  the  “incident  fields”  [48,  49] 

Pine/  \  Vin  o  f  +  kz2t2  COt{kz2z) 

E *  ip'!)  =  77b\s°\ - d ^ - 


£“(/>,  z)  = 


In  ( - 


'  cos \kzJ)  ^ ak J°(bkJ\ 
jkz j€2  +  kz26i  ta n(kz2z) 


(5.44) 


«  =  l  P"Pe 


i(kz2z) 


[Jo(akp)  -  Jo(bkp)) 


(5.45) 


cos(kz2d)  h=\p-pc\ 

where  T)'  and  are  definded  in  Appendix  A. 2,  and  £  is  the  horizontal  distance  between 
the  field  point  and  the  center  of  the  coax  inner  conductor  at  pc.  Once  the  system  (5.33)  is 
solved,  the  input  impedance  of  the  microstrip  patch  antenna  is  found  as 


V 

Z_  VlTl 

in  —  » 


(5.46) 


where  /,„  is  the  current  at  the  base  of  the  coaxial  probe. 

The  computation  of  the  matrix  elements  in  (5.33)  is  extremely  time  consuming,  due 
mainly  to  the  presence  of  Sommerfeld-type  integrals,  which  must  be  repeatedly  evaluated 
for  different  source  and  field  points.  For  the  structure  shown  in  Fig.  5.7,  the  integrals 
encountered  in  Zpp ]  depend  only  on  the  distance  £,  while  those  in  \ZPW ]  and  [ Zwp ] 
depend,  respectively,  on  (£,*')  and  (£,z).  This  suggests  a  time  saving  scheme,  in  which  the 
Sommerfeld  integrals  for  a  given  geometry  and  frequency  are  pre-computed  for  a  grid  of 
points  (using  the  methods  of  Chapter  6),  and  these  values  are  used  in  the  solution  phase  to 
interpolate  the  values  of  the  Sommerfeld  intergals  for  all  required  combinations  of  source  and 
field  points  [50,  17,  18].  We  have  implemented  this  method  using  one- dimensioned  quadratic 
Lagrange  interpolation  for  the  Zpp  terms,  and  two-dimensional  cubic  Lagrange  interpolation 
for  the  Zwp  and  Zpw  terms. 

In  cases  where  the  substrate  is  electrically  thin,  a  simplified  feed  model  may  be  employed 
to  analyze  the  antenna  of  Fig.  5.7,  without  appreciable  loss  of  accuracy.  In  this  model. 
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the  coax  probe  is  replaced  by  a  vertical  filament  of  constant  current  /tn,  which  enters  the 
microstrip  patch  and  spreads  radially  away  from  the  junction  point.  The  excitation  current 
on  the  patch  can  be  approximated  by  IinAJ~(r ),  where  AJ  is  defined  in  (5.30).  This  current 
is  now  considered  to  be  the  source  of  the  “incident  held,”  and  the  current  distribution  on 
the  microstrip  patch  is  the  only  unknown.  Otherwise,  the  solution  procedure  is  similar  to 
that  described  above,  and  it  leads  to  a  matrix  equation 

[ZPP\[IP\  =  [V1*]  (5.47) 


where  Zpp  are  given  in  (5.8),  and  Vp  is  very  similar  in  form  to  ZPJ ,  given  in  (5.34).  Once 
(5.47)  is  solved  for  the  patch  current,  the  antenna’s  input  impedance  is  calculated  as 


Ztn  =  -  f°  Ep  -zdz'  +  Zc 

J  —a 


(5.48) 


where  Ep  is  the  field  due  to  the  total  current  on  the  microstrip  patch  evaluated  at  the 
location  of  the  probe,  and  Zcoax  is  the  correction  term  required  to  account  for  the  “probe 
self-impedance.”  Several  different  expressions  for  Zcoax  have  been  used  in  the  literature  [51, 
52,  53],  with  the  most  elaborate  form  recently  given  by  Mosig  and  Hall  [54,  p.  434]: 


Z^^^-ikody/Tr+j 


.  .  [2d\  a  -  Vo 2  +  4 

irsmh  7  + - u - 


(5.49) 


where  it  is  assumed  that  the  upper  medium  is  free  space  with  the  intrinsic  impedance  tj0, 
and  er  is  the  relative  dielectric  constant  of  the  substrate. 

The  simplified  and  rigorous  coax-feed  models  will  be  referred  to  as  the  1st-  and  2nd-order 
models,  respectively. 


5.4  Transmission  Line  of  Arbitrary  Cross-Section 

The  structure  under  consideration  is  shown  in  a  cross-secliuual  view  in  Fig.  5.9.  The  ground 
plane  at  z  =  —  d  is  assumed  to  be  perfectly  conducting.  The  arbitrarily  shaped  PEC  cylinder 
may  be  open  (for  example,  infinitely  thin  strip)  or  closed  (finite- thickness  strip),  and  is  of 
infinite  extent  and  invariant  along  the  y  coordinate.  As  illustrated  in  Fig.  5.9,  the  contour  of 
the  cylinder  cross-section  is  approximated  in  terms  of  straight  line  segments.  The  definitions 
of  the  ends  of  the  segments  and  the  local  arc-length  coordinate  parallel  to  the  segments 
are  the  same  as  those  in  Section  5.2,  except  that  in  this  case  l  and  r  depend  only  on 
the  transverse  coordinates  x  and  z.  We  assume  that  all  field  components  depend  on  the 
longitudinal  coordinate  y  according  to  e~jk»v,  where  ky  —  0  —  ja  represents  the  propagation 
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Figure  5.9:  Straight  line  segment  model  of  a  cylinder  of  arbitrary  cross-section  embedded  in 
a  grounded  slab. 


58 


constant  of  the  mode.  For  this  two-dimensional  problem,  the  surface  current  density  on  the 
cylinder  can  be  decomposed  as 


J{r)  =  Jt(r)  +  if  Jy(r)  (5.50) 

where  Jt(r)  and  Jy(r )  are  the  transverse  and  longitudinal  surface  current  densities,  respec¬ 
tively,  and  can  be  expressed  as 

J  t[r)  =  Jt{l)e~:kyV  (5.51) 

jy(  T)  =  (5.52) 

where  l  is  the  local  arc- length  coordinate  along  the  circumference  of  the  conductor  and  is 
a  function  of  the  transverse  position  only  (cf.  Fig.  5.9).  Thus,  for  the  transmission  line 
structure,  (5.1)  may  be  expressed  as 


2 


• 

hm  xjr[juAni(l)  +  (Vt  -  yjky)<^m\l)]  =0,  m  =  l,2 

t=i 

(5.53) 

where 

• 

_  •  d  {  .  d 

V,  =  x  —  +  z  — 

OX  oz 

(5.54) 

Am,(l)  =  [  K?A\l\l')-  J{l')dl' 

(5.55) 

L, 

• 

«“(/)  =  /  K?(l\l')q(l')d\'. 

(5.56) 

L, 


The  charge  density  q{l )  in  (5.56)  can  be  obtained  by  using  the  equation  of  continuity, 

?(!)=  X-  (V,  -  yjk,)  ■  J(l).  (5.57) 

The  elements  of  K™  and  K™  can  be  obtained  from  (A.18)-(A.37)  by  making  the  following 


changes: 


So|y(Mj  — *  SoB(ir)]  =  / 

—  oo 

cos  CSi[g{kp)}  — ♦  50(;M(fcr)] 


(5.58) 

(5.59) 


sin C^'i [g{kp)}  — »  ; 


(5.60) 
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£.(r)  = 


,-J^T 


-  k$ 


(5.61 ; 


where  p  =  \J{x  —  x')2  +  (z  —  z')2.  Following  the  MOM  procedure  developed  in  [4]  for  objects 
in  free  space,  the  surface  current  densities  may  be  approximated  as 


N‘ 


j.  =  E£a.(') 


n=  1 


^v  =  E4yn  nU) 


(5.62) 

(5.63) 


n=l 


where  An(l)  is  given  in  (5.17)  and 
'  1,  Ijx— !  <l<ln 

n  „(/)  =  (5.64) 

0,  otherwise  . 

Substituting  (5.55)  and  (5.56),  and  using  (5.62)  and  (5.63)  in  (5.53),  and  testing  the  resulting 
equation  with  An  and  with  y  Iln,  one  obtains  the  matrix  equation 


(5.65) 


7U 
,  pn, 

zty 

pn, 

ft] 

'  0  ‘ 

zyt 

l  prll 

zyy 

.  Pn 

J 

0 

where  the  elements  of  Zpn  are  given  as 
T  h+  . 

=  jw 


r ht  /  ...  .  ....  \  h„ 


4 


z»»  =  ,uht  (a<z++  - 


u> 


in  which 


in 

=  ±  10  (r.n)  =  (‘.y) 


pn  +  ^pr.  y 

(*^‘  -  K) 

(5.66) 

(*^+  -  K+) 

(5.67) 

_  ±. 
K  K 

$>+-) 

P"  1 

(5.68) 

(5.69) 

(t,v)  =  {t,y) 

(5.70) 
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(5.71) 


In 

*in=±  I  K?(i?\r)di'. 


In^  1 

where  ±  in  front  of  the  integrals  (5.70)  and  (5.71)  corresponds  to  the  segment  S*.  In  the 
above  equations,  the  first  and  second  superscripts  ±  correspond,  respectively,  to  segments 
and  S±. 

The  evaluation  of  the  matrix  elements  (c/.(5.70)-(5.71 ))  involves  double  integrations  in 
the  space  and  spectrum  domains.  One  notes  that  the  order  of  the  integrations  may  be 
changed  and  the  integration  in  space  domain  can  be  done  in  closed  form.  Here,  we  only  give 
the  expressions  for  the  matrix  elements  for  the  practically  most  important  case  of  a  strip 
confined  to  the  cover  medium  (region  1): 

<n 


j  |  /')<«'  =  £:  {/“-SJ? 

JnXl 


1 

l 


V\kz2  cot {kz2d)  -  jn2kz\ 


j  K"(lf\l')dl'  =  -~(w~ 
1 

^  n 


inT  1 
In 

J 

fnrfcl 

In 

I 

^n±l 
In 


(5.72) 

(5.73) 

(5.74) 


/  O'? 

^n±  1 

i  n  di1 

(5.75) 

J  KIU‘ V 

ln±  1 

1  n  dr 

(5.76) 

/  +  5“ 

in?l 


eikz2ta.n{kz2d)  +  je2kzl 


Dejkz  i 

+2(/i2^2  —  MlO*^** 


Pn 


D'Dh 


(5.77) 


±± 

pn 


H\kz2  cot(kz2d)  -  jn2kzl 


Dhjk 


+  2(/i2f2  -  A*1  el )‘5’> 


*1 


±± 

pn 


jfczl 


(5.78) 


where  we  have  introduced  the  notations 
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A* 


'“  =  ±7  / 


-  A± 


e~]ktl(zf  +  )e~jkz(xf  -r± ) 


S**  [/(*.)]  =  ±2  J  f(kx) 

—  oo 

sin  [(fc2l  £>**  - 

kzlD^  -  kxD^ 


dkT 


where 


Ppr^  = 


(*f  -  -  Z??ar +  (*?-«- 

D*n±  =  pt-x,  D?=fi-z 


(5.79) 


(5.80) 

(5.81) 

(5.82) 


with  A±  =  /j±/2,  and  where  the  ±  in  front  of  integrals  (5.79)  and  (5.80)  corresponds  to 
the  segment  S*.  Substituting  (5.72)  through  (5.77)  into  (5.70),  and  (5.78)  into  (5.71),  one 
obtains  the  matrix  elements  for  the  case  of  a  strip  confined  to  the  cover  region.  In  a  manner 
similar  to  that  mentioned  in  Section  5.2,  the  integrals  in  (5.70)  and  (5.71)  may  be  reused  for 
different  matrix  elements. 

The  homogeneous  system  (5.65)  has  non-trivial  solutions  for  those  values  of  ky,  which 
render  its  determinant  vanish.  Hence,  to  obtain  the  propagation  constants  of  the  various 
modes  of  the  microstrip,  a  search  is  performed  for  the  zeroes  of  the  determinant  in  the 
complex  fcy-plane.  For  each  propagation  constant,  the  homogeneous  system  (5.65)  is  solved 
for  the  corresponding  modal  current  distribution. 
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Chapter  6 

Numerical  Evaluation  of  Sommerfeld 
Integrals 


The  dyadic  and  scalar  kernels  of  the  MPIE  comprise  Sommerfeld  integrals  Sn(f)  of  the  form 
(3.3).  When  the  MPIE  is  solved  by  the  MOM,  as  discussed  in  the  previous  chapter,  these 
integrals  must  be  numerically  evaluated  for  wide  ranges  of  variation  in  media  parameters, 
frequency,  and  spatial  variables.  Several  numerical  algorithms  for  the  computation  of  those 
integrals  are  developed  in  this  chapter,  based  on  the  deformation  of  the  integration  path  off 
the  real  axis  into  the  complex  -plane.  Section  6.1  contains  a  detailed  discussion  of  the 
integrals  associated  with  the  half-space  problem.  One  of  the  techniques  developed  therein  is 
also  employed  to  evaluate  the  improper  integrals  that  arise  in  the  grounded  slab  geometry 
pertinent  to  microstrip  patch  antennas.  Section  6.2  is  devoted  to  the  evaluation  of  spectral 
integrals  associated  with  the  microstrip  transmission  line  problem. 

6.1  Sommerfeld  Integrals  for  the  Half-Space  Problem 

The  Sommerfeld  integrals  encountered  in  this  work  can  be  expressed  as  ( cf .  (3.3)) 

OO 

$.[/(*,)]«  /  f{k,)  J,((  k,)k,'*‘ dk,.  (6.1) 


For  the  half-space  problem,  the  integrand  in  (6.1)  is  given  as 

7/i  )  N(M,k2;kp)  f  e~jfc,l|z+z'1,  m  =  i 

I{  p)  ~  D(kuk2\kp)  {  m^i 


(6.2) 


where  the  functions  N  and  D  can  assume  different  forms  (cf.  Appendix  A.l),  and  where 

(6.3) 


kti  =  \Jk?  -  kp2,  i  =  1,2. 
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The  function  kz{  has  branch  points  at  kp  =  ±kt.  For  analytical  convenience,  we  will  assume 
that  the  media  are  lossy;  that  is,  kt  =  k[  —  jk",  where  both  k[  and  k "  are  positive  numbers.  A 
lossless  medium  will  be  considered  to  be  the  limit  of  a  lossy  medium  as  the  dissipation  goes 
to  zero  ( k ”  — ►  0).  We  also  assume  that  k'2  >  k[  and  k2  >  k ",  which  hold  in  most  practical 
situations.  To  specify  kzl  uniquely,  we  may  view  the  complex  kp- plane,  where  kp  =  ('  + 
as  a  two-sheeted  Riemann  surface  with  the  sheets  connected  along  the  branch  cuts.  On 
this  surface,  kzi  is  a  single- valued  analytic  function  of  kp.  The  choice  of  the  branch  cuts  in 
the  complex  ^-plane  is  rather  arbitrary,  although  a  particular  choice  may  be  convenient 
for  a  specific  problem.  In  this  work,  we  exclusively  use  the  branch  cuts  specified  by  the 
requiremert  that  Im(kzt)  <  0  everywhere  on  one  sheet,  referred  to  as  the  top  or  proper 
sheet.  The  other  sheet,  on  which  Im (kzt)  >  0,  is  then  referred  to  as  the  bottom,  or  improper 
sheet.  This  definition  implies  that  the  two  sheets  are  joined  together  by  the  curve  defined 
by  Im(fczt)  =  0,  which  therefore  locates  the  desired  branch  cut.  Let  us  first  write  [55] 

kzi  =  |  kzx  |  *eJ*  =  k2  -  kp 2,  -  2tt  <  9  <  2n  (6.4) 

whence 

fczi  =  |  kzt  |  ej9/2.  (6.5) 

It  is  clear  that  if  Im(A;ZI)  <0,-7 r  <  9/2  <  0,  and  if  lm(kz,)  >  0,  0  <  9/ 2  <  7r.  This  suggests 
that  we  define  a  two-sheeted  A:^-plane,  where  —  2ir  <  9  <  0  on  the  top  sheet  and  0  <  9  <  2x 
on  the  bottom  sheet,  as  shown  in  Fig.  6.1.  These  sheets  are  joined  along  the  positive  real 
axis,  where 

Im(fc^)  =  0  and  Re(^)  >  0.  (6.6) 

For  later  reference,  we  also  show  in  Fig.  6.1, the  regions  where  Re(fcz,)  >  0  and  R e(kzi)  <  0  in 
the  entire  two-sheeted  kzi- plane.  We  now  map  the  two-sheeted  k2t -plane  into  a  two-sheeted 
kp2- plane  by  means  of  (6.4).  The  result  is  shown  in  Fig.  6.2,  where  we  have  introduced  the 
notation 

k2  =  (k?  -  k'/2)  -  j2k\k”  =  r  +  ;U  (6.7) 

It  follows  from  (6.4)  and  (6.7)  that  the  branch  cuts  in  the  kp2- plane  are  defined  by 

Re(fcp2)  <  t  and  Im(fcp2)  =  ft.  (6.8) 
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R e(fc„)  <  0  R e(kzl)  <  0 

bottom  sheet  (Im (ktl)  >  0) 

(b) 


Figure  6.1:  Two  sheets  of  the  k^- plane,  (a)  Top  sheet  (Im£n  <  0)  and  (b)  bottom  sheet 
(Im kzt  >  0). 
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Figure  6.2:  Two  sheets  of  the  fc^-plane.  (a)  Top  sheet  (ImA:z,  <  0)  and  (b)  bottom  sheet 
(Im kz,  >  0). 
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Finally,  the  two-sheeted  complex  fcp2-plane  is  mapped  into  a  two-sheeted  kp- plane.  Noting 
that  kp 2  =  C2  -  ("2  +  j2('(",  we  conclude  from  (6.8)  that  the  branch  cuts  in  the  fcp-plane 
are  specifed  by 

C'2  -  C"2  <  r,  and  2('C  =  V.  (6.9) 

This  mapping  is  illustrated  in  Fig.  6.3. 

So  far,  we  have  only  considered  a  single  square  root  function  kzi.  However,  the  integrand 
of  (6.1)  depends  on  both  kz t  and  kz 2.  Therefore,  to  make  it  single- valued,  we  must  introduce 
two  pairs  of  branch  cuts  in  the  fcp-plane,  defined  by  Im {kzi)  =  0  and  Im (kz2)  =  0.  respec¬ 
tively  (Fig.  6.4a).  We  may  now  view  the  kp- plane  as  a  four-sheeted  Riemann  surface,  as 
illustrated  schematically  in  Fig.  6.46  [33].  Observe  that  everywhere  on  Sheet  I,  Im(/bzi)  <  0 
and  Im(A;Z2)  <  0. 

In  addition  to  the  branch  point  discussed  above,  the  integrand  of  (6.1)  may  exhibit  a 
pair  of  first  order  poles  at  the  zeroes  of  the  denominator  function 


D{k\,  k2\  kp)  =  t\kZ2  +  e2kzl 


(6.10) 


which  are  given  by 


kpp  —  di  , 


k\k\ 

i|*rr*r 


(6.11) 


These  poles  only  exist  on  Sheets  I  and  IV.  Their  location  is  indicated  in  Fig.  6.4a,  which 
represents  Sheet  I  of  the  Riemann  surface. 

In  view  of  the  above  discussion  it  should  be  clear  that  to  ensure  the  convergence  of  the 
integral  (6.1),  the  integration  path  should  be  selected  on  Sheet  I,  as  illustrated  in  Fig.  6.4a. 

Numerical  evaluation  of  the  Sommerfeld  integrals  (6.1)  is  difficult  because  of  the  oscilla¬ 
tory  behavior  of  the  integrand  and  its  rapid  variation  near  the  singularitites  (branch  points 
and  poles).  Various  numerical  integration  techniques  have  been  developed  to  carry  out  these 
integrals.  The  real-axis  path  has  been  used  in  related  problems  by  Siegel  and  King  [56],  Kuo 
and  Mei  [57],  Lin  and  Mei  [58],  Katehi  and  Alexopoulous  [59],  Johnson  and  Dudley  [60], 
and  Mosig  and  Gardiol  [17,  61].  The  techniques  of  deforming  the  integration  path  to  a  con¬ 
tour  off  the  real  axis  to  avoid  the  singularities  and  to  accelerate  convergence  of  the  integrals 
have  been  developed  by  Miller  et  al.  [62],  Sarkar  [63],  Burke  et  al.  [50],  Michalski  [24]  and 
Michalski  and  Butler  [64].  The  method  of  deforming  the  integration  path  from  the  real  axis 
to  vertical  branch  cuts  has  been  used  by  Fuller  and  Wait  [65],  and  Kong  et  al.  [66].  Several 
approximate  approaches  have  been  used  by  Hansen  [67],  and  by  Mosig  and  Gardiol  [61]. 
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Figure  6.3:  Two  sheeted  fcp-plane.  (a)  Top  sheet  (Im <  0)  and  (b)  bottom  sheet  (Imfcr,  > 

Figure  6.4:  (a)  The  &p-plane  showing  the 
surface  of  four  sheets  in  the  A^-plane. 


Im(^i)  <  0,  hn(kx2)  <  0 
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anch  cuts  for  kz i  and  (b)  The  Riemann 


The  steepest  descent  path  (SDP)  is  probably  the  most  efficient  one  to  use  [24],  However, 
since  a  closed  form  expression  for  the  SDP  is  not  available  when  r  and  r'  are  in  different 
layers,  we  employ  in  this  work,  as  a  compromise,  four  suboptimal  paths,  which  are  chosen 
depending  on  the  relative  values  of  =  |  p  —  p'  f  and  h  =  |  z  |  +  |  z'  | ,  and  on  the  param¬ 
eters  of  the  media.  The  necessity  to  evaluate  Bessel  functions  of  complex  arguments  on 
these  paths  does  not  significantly  affect  the  efficiency  of  our  approach,  because  it  has  been 
demonstrated  by  Grun  and  Rahmat-Samii  [68]  that  the  simple  polynomial  approximations 
of  Bessel  functions  given  in  [69]  can  be  analytically  continued  into  the  complex  plane.  Since 
all  the  paths  employed,  or  similar,  have  been  previously  described  in  the  literature,  we  only 
summarize  them  below. 

Path  I 

The  integration  path  developed  by  Michalski  and  Butler  [64]  is  illustrated  in  Fig.  6.5a. 
This  path  avoids  the  real-axis  singularities  and  is  the  SDP  for  the  exponential  function  part 
of  the  integrand  (but  not  for  the  whole  integrand).  Hence,  on  this  path  the  exponential 
function  decreases  monotonically,  because  the  imaginary  part  of  the  exponent  is  by  definition 
constant  on  the  SDP.  When  the  source  and  observation  points  are  in  the  same  layer  i  for 
i  =  1,2,  one  can  use  the  transformation 


kp  —  sk0yj2jnt  +  s2,  0  <  s  <  oo 

(6.12) 

where  s  is  a  real  variable,  n,  =  and  Re^/-  >  0.  If  i  ^  m,  i.e 

are  in  different  layers,  we  choose 

.,  the  source  and  field  points 

/  -2  C 

(6.13) 

k0  ~  koK  - -======,  0  <  s  <  oo 

p  V  B  +  y/B2-4AC 

in  which  the  ReyC  >  0  and  Imy^”  <  0  branches  should  be  selected  for  the  external  and 

internal  roots,  respectively.  The  parameters  in  (6.13)  are  given 

as 

A  =  (z'2  -  z2)2 

(6.14) 

B  =  4|  z' z  j  de  —  2s(z'2  +  z2) 

(6.15) 

C  =  .«2  —  4sn,nm|  z  z  \ 

(6.16) 

with 


d  =  n,|  z  |  +  nm  |  z  | 


(6.17) 
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R e(fcp) 


Re(fcp) 


e  =  n,|  2  |  +  nm|  z 


(6.18) 


Since  the  oscillations  due  to  the  Bessel  function  are  still  present,  this  path  selection  is  most 
efficient  when  h  =  |  a'  |  +  |  ?  |  >£.  Also,  the  exponential  growth  of  the  Bessel  function  away 
from  the  real  axis  limits  how  far  the  path  can  veer  off  into  the  complex  plane,  which  restricts 
the  applicability  of  this  method.  In  ffie  case  m  =  i.  this  method  is  applicable  provided  £  is 
electrically  sma'1  in  the  mediumi,  be.. 


1 


(6.19) 


When  m  ^  z,  this  restriction  becomes 
< 


1  z'  \ 
1  z  1 

+ 

1 

n(| 

l*'l 

+  n'  1 

1  m ! 

z 

(6.20) 


To  numerically  evaluate  the  integral  i)6.1)  along  this  path,  a  low-order  Gauss- Legrandre 
quadrature  is  repeatedly  employed  between  consecutive  “break  points,”  which  are  derived 
from  the  real-axis  zeroes  of  the  Bessel  function.  When  each  of  three  consecutive  intervals 
contributes  less  than  10~4  of  the  accumulated  sum,  the  semi-infinite  integration  is  terminated 
end  is  said  to  converge. 

Path  II 

An  alternate  path  [50],  which  initially  sprays  from  the  real  axis  to  avoid  the  branch  points 
at  and  k2,  then  continues  parallel  to  the  real  axis,  as  shown  in  Fig.  6.5 6,  has  been  used 
to  overcome  the  limitations  (6.19)  and  (6.20)  in  using  Path  I.  The  parameter  P  is  chosen 
according  to 


P  =  minQ,i).  (6.21) 

which  restrets  the  contour  to  small  values  of  £  lm(kp),  to  limit  the  growth  of  the  Bessel 
function.  This  compromise  path  accelerates  the  convergence  of  the  exponential  function 
(even  though  it  is  not  the  SDP),  while  avoiding  the  overflow  in  the  evaluation  of  the  Bessel 
function.  On  this  path  we  employ  the  same  numeibal  integration  technique  as  that  used  on 
Path  1 

Path  III 

Whx'u  (  >  /i.  especially  when  h  =  0,  Paths  I  and  II  are  net  efficient,  because  there  may 
be  many  oscillations  of  he  Bessel  function  before  convergence  is  achieved.  Therefore,  we 
introduce  Path  III.  which  in;tia!ly  strays  from  the  real  axis  to  avoid  the  singularities,  then 
returns  at  A  .  =  T  to  th  .  real  axis.  The  point  T  is  arbitrarily  chosen  as  T  =  A .> ( ri2  +  1) 


i  .i 


(recall  that  we  assume  k'2  >  &().  On  the  first  part  we  employ  the  same  numerical  integration 
technique  as  that  used  on  Path  I.  On  the  real  axis  (kp  >  T)  we  employ  the  method  of 
averages  [17,  61,  70],  which  is  summarized  below. 

Let  us  consider  the  integral 


OO 

/  =  J  f(kp)JnUkp)kpn+ldkp 


(6.22) 


where  f{kp)  is  a  continuous  function,  and  Jn(£kp)  has  successive  zeroes,  kPm ,  superior  to  the 
integration  boundary  T.  It  follows  from  the  method  of  averages  [17,  61]  that 


/%/,M  =  21~mT  (M  M/I 


m= 1 


m  —  1 


(6.23) 


where 


nPm 

il  =  J  f{kP)  Jn(£kp)kpn+1  dkp,  m=  1,2,..., 


M. 


(6.24) 


The  choice  of  M  in  (6.23)  depends  on  the  requirement,  of  the  accuracy  in  the  problem.  In  this 
work,  we  choose  M  =  5.  The  method  of  averages  is  especially  suitable  when  the  observation 
and  source  points  are  on  the  interface.  We  find  that  one  still  can  obtain  the  accurate  result 
of  Sommerfeld  intcgra's  by  using  this  method  when  the  observation  and  source  points  are 
off  the  interface,  provided  that  £  h. 

Path  IV 

integration  along  Path  III  becomes  inefficient  when  <f  becomes  large,  bacause  there  may 
be  many  oscillations  of  the  Bessel  function  between  0  and  T .  In  this  case,  it  is  useful  to 
deform  the  integration  path  to  two  vertical  paths,  as  shown  in  Fig.  6.5d.  For  this  purpose, 
by  using  the  relation 

2Jn(x)  =  H^{x)  +  (-1  r+lHl2)(-x)  (6.25) 

where  H^'  is  the  Hankel  function  of  the  second  kind  and  order  n,  one  can  express  (6.1)  as 

Sn ( f)  =  ~  j  f(kzl,kl2-kp)Qn^kp)Kpn+le~^dkp  (6.26) 

—  OO 

where  we  ^ave  '  troduced 

Qn(ik„ )  =  H™tikp)e*k>  (6.27) 
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which  has  the  asymptotic  form  [69] 


(6.28) 


It  is  noted  the  transformation  (6.25)  introduces  a  branch  point  singularity  at  kp  =  0  and  the 
associated  branch  cut  on  the  negative  real  axis  (Fig.  6.5d).  Hence,  the  integral  (6.26)  should 
be  taken  in  the  Cauchy  principal  value  sense  for  n  >  1.  As  a  result,  (6.26)  becomes 


Sn(f)  =  J g{kp)  dkp  -;7rRes{tf(0)} 


(6.29) 


where 


9(W  =  ^J(h)QMkr)k,’'+'e-*'” 


(6.30) 


and  Res  {#(0)}  is  the  residue  of  the  integrand  at  kp  =  0.  The  integration  path  C  is  that 
illustrated  in  Fig.  6.5d.  In  the  following  analysis,  we  will  refer  to  this  path  as  the  “original” 
path  of  integration. 

We  observe  that  when  Im(fc„)  becomes  negative,  the  integrand  in  (6.26)  decays  exponen¬ 
tially.  It  is  clear  that,  for  large  values  of  £,  one  should  deform  the  integration  path  into  the 
vertical  paths  [65,  66]  as  illustrated  in  Fig.  6.5d.  The  integral  (6.26)  can  be  written  in  this 


case  as 


sn  =  Si1*  +  S<2> 


(6.31) 


where  S^1*  and  S^2*  are  the  integrals  along  the  two  vertical  paths  emanating  from  k\  and 
k2,  respectively,  as  shown  in  Fig.  6.5d.  Referring  to  Figs.  6.3  and  6.4,  we  observe  that  the 
integration  path  C i  starts  on  Sheet  IV,  where  Re(kz\,kz2)  >  0,  and  proceeds  upward  along 
the  vertical;  it  leaves  Sheet  IV  when  it  crosses  the  kz2  branch  cut,  and  enters  a  region  of 
Sheet  II,  where  Re(kzX, kz2)  >  0;  it  then  reaches  Sheet  I,  where  R e(fc2l)  <  0  and  R e(kz2)  >  0, 
when  it  crosses  the  kzi  branch  cut;  it  turns  around  the  branch  point  at  k\  and  goes  down 
along  the  vertical;  finally,  it  leaves  Sheet  I  when  it  crosses  the  kz2  branch  cut,  and  enters  a 
region  of  Sheet  III,  where  Re(kzi)  <  0  and  Re(A:r2)  >  0.  Similar  analysis  indicates  that  the 
path  C2  starts  on  Sheet  III,  proceeds  upward  along  the  vertical  from  Sheet  III  to  Sheet  I, 
and  reaches  the  branch  point  at  k2.  On  this  part  of  C2 ,  Re(kzi)  <  0  and  Re(fcr2)  0-  The 
integration  path  C2  finally  emanates  from  k2  and  proceeds  down  along  the  vertical  on  Sheet  I, 
where  Re(kzl,  kz2)  <  0.  From  the  above  considerations,  one  can  determine  the  branches  of 
the  square  root  functions  kzl  and  kz2  along  the  vertical  paths  as  follows: 
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(6.32) 


Rc(  kz\  ,  kz2^  0,  on  Ci 

Re(kz\)  <  0,  Re(£j2)  >  0,  on  Cf  and  C2~ 
Re(/tzi ,  kz2)  <  0,  on  C% . 

Introducing  the  variable  transformations 
kp  —  k\  —  jk0s2,  0  <  s  <  00,  on  Ci 
kp  =  k2  —  jkos2,  0<s<  00,  on  C2 
we  obtain  the  following  integrals  by  invoking  (6.32): 

S<x)  =  -—  J  [f(-kzi,kz2;  kp)  -  f{kzUkz2\ kp) 


(6.33) 

(6.34) 


■Qn(ikp)k;+1 


(-2 jk0s)e~koi32  ds 
kp  =  k0(ni  -  js 2) 


(6.35) 


and 


•<?.«*,)  C1 


(-2;3t05)e-^^2  ds 
'fcP  =  ^0(n2  -  js2) 


(6.36) 


where  it  is  assumed  that  Re(kzl,  kz2)  >  0.  The  proximity  of  the  pole  at  kpp  (cf.  (6.11)) 
on  Sheet  I  to  the  integration  path  Ci  makes  the  integrand  vary  rapidly  near  s  =  0.  To 
remedy  this  difficulty,  a  pole  subtraction  technique  is  employed  [71],  which  transforms  the 
integral  (6.35)  as 

|  /  f  r  ~  ~  i 

J  \f{-kzi,kz2]kp)  -  f{kzUkz2-,kp)\ 


si1*  = 


■QMkP)kp+1 


-2 jkos)  -  e"^'3  ds  +  jTrRpW  (0^3,)  J  (6.37) 

kp  =  k0(ni  -  js2) 

where  Ftp  is  the  Residue  of  {f{kzi,  kz2\ kp)  Qn(£kp)kp+l}  at  kp  =  kpp,  sp  is  the  pole  location 
in  the  s-plane,  and  W(z)  =  e“*3erfc(  —  jz)  is  the  Error  function  [69]. 

Based  on  the  considerations  above  and  on  our  computational  experience,  a  set  of  guide¬ 
lines  has  been  developed  for  determining  the  “best”  path  of  integration.  The  guide  chart  is 
given  in  Fig.  6.6. 
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The  integrals  that  arise  for  the  grounded  slab  geometry  pertinent  to  microstrip  patch 
antennas  are  similar  to  those  considered  above.  Their  integrand  functions  have  similar 
behavior  to  those  of  the  half-space  problem,  except  that  they  do  not  have  branch  points 
associated  with  the  slab  region,  but  have  surface  wave  poles,  which  are  located  betwen  k\ 
and  in  the  fcp-plane.  To  evaluate  those  integrals,  we  have  employed  Path  III  shown  in 
Fig.  6.5c,  in  conjunction  with  the  interpolation  scheme  described  in  Section  5.3. 


6.2  Sommerfeld  Integrals  for  the  Transmission  Line 
Problem 


In  this  section,  to  simplify  the  discussion,  we  only  consider  the  case  where  the  infinite 
conductor  is  confined  to  the  cover  medium  (region  1).  A  typical  integral  that  arises  in  this 
context  is  given  in  (5.80).  This  integral  can  be  expressed  in  terms  of  the  integrals  of  the 
form 


j  _  7  N(k*uk*ikx)  ikwt(t+tl)  \  cos[fc*(x  -  x')]  1  , 

J  DeDh  1  sin[£x(x  -  x')]  J 


(6.38) 


where  the  cosine  or  sine  functions  arise  depending  on  whether  the  numerator  functions  N 
are  even  or  odd  in  kx.  The  denominator  functions  De  and  Dh  are  given  in  (A. 38)  and  (A. 39), 
respectively,  and 


Ki  =  \A?  ~kl~  i  =  1 ’2- 


(6.39) 


The  propagation  constant  will  be  expressed  as 


ky  =  0-ja  (6.40) 

where  0  is  the  phase  constant  and  a  is  the  attenuation  constant.  It  will  be  assumed  that 
a  >  0  and  0  >  0. 

It  can  be  shown  that  the  integrand  of  (6.38)  is  an  even  function  of  kzi.  Therefore, 
unlike  in  the  half-space  problem,  only  one  pair  of  branch  points,  associated  with  kzl ,  arise 
in  the  present  case.  These  branch  points  and  the  associated  branch  cuts,  selected  by  the 
criteria  given  in  Section  6.1,  define  a  fcr-plane  as  a  two-sheeted  Riemann  surface.  As  before, 
Im(£zi)  <  0  on  the  top  sheet  and  the  opposite  holds  on  the  bottom  sheet.  Since  the  value 
of  ky  will  be  specified  on  each  iteration  of  the  search  procedure,  we  express  kz\  as 

kzi  =  \Jk2  ~  kl  (6-41) 
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where  we  have  defined 


—  k2  —  k2  —  (&2  —  02  +  a2 )  4-  j2a0 

—  t  +  ?H. 


(6.42) 


We  assume,  for  simplicity,  that  the  media  are  lossless.  Consequently,  ki  in  (6.42)  is  real. 

With  the  notation  introduced  above,  we  can  directly  use  the  analysis  which  led  to  Fig.  6.2, 
provided  that  A^2-plane  is  replaced  by  k2-plane.  Therefore,  referring  to  Fig.  6.2,  we  define 
the  branch  cut  in  the  k2- plane  by 


Re(A:2)  <  r  and  Im(A;2)  =  0. 


(6.43) 


These  relations  define  a  pair  of  branch  cuts  in  the  A^-plane,  whose  shape  depends  on  the 
values  of  ki  (and  thus  frequency)  and  ky  (or  a  and  0).  In  addition  to  the  branch  points, 
there  will  also  be  a  finite  number  of  poles  in  the  top  sheet  of  the  itr-plane,  contributed  by 
the  zeroes  of  the  denominator  functions  De  and  Dh .  These  poles  are  easily  determined  in  the 
A;, -plane  (recall  that  kp 2  =  k2x  +  A;2),  where  they  are  located  on  the  real  axis  between  k\  and 
k2,  and  also  between  ~k2  and  —kx  [32].  Their  locations  in  the  Arr-plane  will  depend  on  the 
value  of  ky.  In  the  following  analysis,  we  assume  for  simplicity  that  only  one  pair  of  poles, 
±ks,  associated  with  the  TM0  mode  of  the  slab,  occurs.  This  is  the  case  most  frequently 
encountered  in  practice. 

As  mentioned  above,  the  shape  and  location  of  the  branch-cuts  in  the  A:x-plane  depends 
on  the  value  of  ky.  Following  the  analyses  of  [72,  73,  74],  we  divide  the  range  of  ky  into  three 
subregions,  according  to 


(f)  k2  <  02  <  k2,  a  =  0 
(«)  k2  <-0 2  -  q2  <  k] 
(ni)  — oo  <  02  —  q2  <  k\ 


(6.44) 


where  02  —  a2  =  Re(fc2).  In  region  i,  ky  is  real  and  the  stripline  mode  is  bound,  t'.e.,  it 
propagates  unattenuated.  The  fields  of  bound  modes  are  concentrated  in  the  region  of  the 
strip  and  decay  exponentially  in  transverse  directions  away  from  it.  In  regions  ii  and  in',  kv 
is  complex  (a  >  0),  and  the  mode  is  said  to  be  leaky  [72,  73,  74].  In  this  regime,  the  mode  is 
attenuated  as  it  propagates,  due  to  the  leakage  of  energy  into  the  environment.  In  region  ti, 
there  is  only  leakage  into  a  surface  wave  of  the  slab,  which  propagates  away  from  the  strip. 
It  can  be  shown  [72]  that  the  field  of  such  mode  decays  exponentially  in  the  z  direction, 
but  increases  exponentially  away  from  the  strip  in  ±x  directions.  For  ky  in  region  tit,  the 
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leakage  is  both  into  the  surface  wave  and  into  the  space  wave.  The  field  amplitude  increases 
exponentially  in  the  directions  perpendicular  to  the  strip  [72,  75]. 

Let  us  first  consider  the  bound  modes,  for  which  ky  is  in  region  i  (c/.(6.44)).  With  r  <  0 
and  Q  =  0  in  (6.42),  the  two-sheeted  complex  A^-plane  is  mapped  into  a  two-sheeted  kx- 
plane  as  illustrated  in  Fig.  6.7a.  The  poles  are  in  this  case  located  in  the  ( —  i  J~k^  —  k\. 
j\Jky  —  k\)  interval  of  th.2  imaginary  &r-axis,  as  shown  in  Fig.  6.7a.  The  integration  path  in 
(6.38)  is  chosen  along  the  real  axis  on  the  top  sheet  (Fig.  6.7a). 

When  ky  passes  from  region  i  to  region  it,  ky  acquires  the  imaginary  part,  so  that  r  <  0 
and  0.  >  0.  In  this  process  the  pole  on  the  negative  imaginary  axis  in  the  fcr-plane  moves 
to  the  origin,  then  enters  the  first  quadrant,  as  illustrated  in  Fig.  6.7 6.  In  a  like  manner, 
the  pole  on  the  positive  imaginary  axis  moves  to  the  third  quadrant.  The  integration  path 
must  be  deformed  so  that  the  pole  in  the  first  quadrant  lies  below  the  path,  as  illustrated  in 
Fig.  6.76.  The  first  quadrant  pole  contributes  an  exponentially  increasing  surface  wave,  in 
agreement  with  the  discussion  above. 

When  ky  moves  from  region  ii  to  region  iii,  the  value  of  r  changes  sign,  hence  in  region  iii 
r  >  0  and  f!  >  0  (cf.  (6.42)  and  (6.44)).  At  the  transition  point  r  =  0,  the  branch  points 
in  the  first  and  third  quadrants  switch  positions  along  the  diagonal  line  in  the  A:r-plane. 
Therefore,  the  integration  path  must  go  above  the  branch  point  k  in  the  first  quadrant,  as 
illustrated  in  Fig.  6.7c.  One  can  prove  that  this,  in  turn,  leads  to  a  modal  field  which  grows 
exponentially  in  the  z  direction  [75].  The  dashed  line  on  the  integration  path  in  Fig.  6.7 c 
signifies  that  the  path  is  on  the  bottom  sheet. 

As  the  above  analysis  indicates,  different  integration  paths  must  be  employed  depending 
on  which  mod'1  is  being  computed.  The  various  paths  are  suggested  in  Fig.  6.7.  In  the 
computer  implementation,  to  accelerate  the  convergence  of  the  integral  and  to  avoid  the 
singularities,  slightly  modified  paths,  shown  in  Fig.  6.8,  have  been  employed.  For  bound 
modes,  we  use  the  paths  shown  in  Fig.  6.8a,  where  the  value  of  P  is  arbitrarily  chosen  as 
P  —  k\.  For  ky  in  region  ii,  the  paths  in  Fig.  6.86  are  employed,  where  P  is  the  same  as  that 
in  Fig.  6.8a  and  the  value  of  T  is  chosen  as  T  =  \  Re(fcf  —  k*)%  |  +  ki,  to  ensure  that  the 
pole  is  located  to  the  left  of  this  point.  The  integration  around  the  pole  is  carried  out  by 
the  method  of  residues.  For  ky  in  region  iii,  we  employ  the  paths  shown  in  Fig.  6.8c,  where 
T  is  the  same  as  that  in  Fig.  6.86  and  the  value  of  P  is  chosen  as  P  =  |  Im(A:j  —  k %)?  |  +  k\, 
to  ensure  that  the  branch  point  in  the  first  quadrant  lies  below  the  paths.  In  all  three  cases, 
when  z  =  z'  =  0  (this  is  always  the  case  for  an  infinitely  thin  strip),  we  use  path  C\  in 
Fig.  6.8.  When  the  path  is  on  the  real  axis,  the  method  of  averages  [17]  discussed  in  the 
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previous  section  is  employed  to  accelerate  the  rate  of  convergence.  When  h  =  |  z  \  +  ]  z'  |  >  0 
path  C2  is  used  to  benefit  the  convergence  from  the  exponential  function  associated  with  h 


83 


Chapter  7 

Numerical  Results 


In  this  chapter,  we  present  sample  numerical  results  for  antennas  and  scatterers  partially 
buried  in  earth  or  water,  for  coax-fed  microstrip  patch  antennas,  and  for  open  microstrip 
lines.  In  all  cases,  the  medium  of  region  1  is  free  space  with  parameters  ej  =  e0  and  /zi  =  /i0. 
The  medium  of  region  2  is  characterized  by  e 2  and  /i2  =  fi0,  and  it  may  be  lossy  or  lossless. 
In  the  former  case,  its  permittivity  is  complex  and  is  given  as  e2  =  fo^r  —  j&/u,  where  er 
and  a  are  the  relative  dielectric  constant  and  the  conductivity  of  the  medium,  respectively, 
and  u)  =  27r/.  The  loss  tangent,  tan  6  =  cr/(ue),  is  also  used  to  characterize  the  losses  of  the 
medium. 

In  Sections  7.1  and  7.2,  to  facilitate  the  interpretation  of  the  results,  we  take  in  all 
cases  /  =  300  MHz,  which  corresponds  to  free  space  wavelength  A0  =  1  m.  The  results  for 
driven  antennas  assume  unit-strength  <5-gap  generators,  and  the  results  foi  scanerers  assume 
illumination  by  plane  waves  incident  normally  on  the  interface. 

7.1  Surfaces 

We  first  consider  a  straight,  inclined,  thin-wire  antenna  with  radius  a,  partially  buried  in  dry 
earth,  as  illustrated  in  Fig.  7.1a.  In  the  numerical  procedure,  the  wire  was  approximated 
by  a  flat,  narrow  strip  of  equivalent  width  4a  [76],  which  was  modeled  by  60  triangular 
patches,  as  shown  in  Fig.  7.16.  In  Fig.  7.2,  the  computed  current  distributions  on  the  wire 
are  compared  with  the  data  from  the  NEC  [27].  A  good  agreement  is  observed  for  the 
inclination  angle  a  =  45°  (Fig.  7.2 a),  while  for  a  =  80°  (Fig.  7.26)  the  agreement  is  less 
favorable.  A  possible  explanation  of  this  slight  disagreement  is  the  fact  that  in  the  NEC 
formulation  a  charge  discontinuity  condition  at  the  interface  is  enforced,  which  is  strictly 
valid  only  for  the  vertical  antenna  [27].  Obviously,  this  condition  would  affect  the  solution 
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Figure  7.1:  (a)  Inclined  thin-wire  antenna  partially  buried  in  earth  and  (b)  its  strip  model 
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more  for  a  =  80°  than  for  a  =  45° . 

In  Fig.  7.3  is  shown  a  relatively  thick,  vertical,  cylindrical  antenna,  which  penetrates  the 
interface  between  two  media.  This  so-called  ground  stake  antenna  was  previously  analyzed 
by  Butler  and  Michalski  [77].  In  the  numerical  procedure,  the  circular  cylinder  was  approxi¬ 
mated  by  a  square  cylinder  with  the  same  circumference,  and  was  modeled  by  112  triangular 
patches.  We  show  the  axial  current  distribution  on  the  antenna  for  the  cases  where  the  lower 
medium  is  dry  earth  (Fig.  7.4a)  and  salt  water  (Fig.  7.46).  Our  results  are  seen  to  perfectly 
agree  with  the  data  obtained  in  [77].  Since  our  code  can  handle  surfaces  of  arbitrary  shape, 
it  was  a  simple  matter  to  examine  the  effect  of  putting  the  end  caps  on  the  hollow  cylinder. 
The  resulting  current  distributions  are  also  displayed  in  Fig.  7.4.  One  observes  that,  as 
expected,  the  caps  only  affect  the  current  near  the  ends  of  the  cylindrical  surface. 

The  results  in  Fig.  7.5  are  for  a  finite,  hollow,  horizontal  cylinder,  which  is  partially  buried 
in  a  dielectric  medium  (see  the  inset).  The  corresponding  two-dimensional  problem  has  been 
solved  by  Xu  [26],  who  employed  the  magnetic  field  integral  equation  formulation.  One 
observes  that  both  the  magnitude  (Fig.  7.5a)  and  phase  (Fig.  7.5 6)  of  the  normalized  current 
distribution  along  the  circumference  and  in  the  middle  of  the  finite  tube  agree  favorably  with 
the  corresponding  results  for  the  infinite  cylinder.  In  the  triangle-patch  model  of  the  cylinder, 
288  patches  were  employed. 

In  Figs.  7. 6-7.9,  we  present  results  for  a  flat,  rectangular  plate  partially  buried  in  dry  earth 
(see  the  inset  in  Fig.  7.6a).  The  magnitudes  of  the  dominant  and  transverse  components  of 
the  current  distribution  on  the  plate  are  shown  in  Fig.  7.7 a  and  Fig.  7.76  for  the  inclination 
angle  <y  =  30°,  and  in  Fig.  7.9a  and  Fig.  7.96  for  a  =  60°.  In  the  triangle-patch  model  of  the 
plate,  252  patches  were  employed. 

The  slight  nonsymmetry  with  respect  to  the  center  line  of  the  plate  observed  in  the 
transverse  component  of  the  current  is  due  to  the  fact  that  the  symmetry  of  the  plate  was 
not  preserved  in  its  triangle- patch  model.  In  Figs.  7.6  (a  =  30°)  and  7.8  (a  =  60°)  we 
compare  the  dominant  component  of  the  current  distribution  along  the  center  line  of  the 
plate  with  the  corresponding  result  for  an  infinite  strip  [26].  One  observes  a  reasonably  good 
agreement  between  the  two  results,  both  in  magnitude  (Fig.  7.6a)  and  phase  (Fig.  7.66)  for 
a  =  30°,  and  a  less  favorable  agreement  for  a  =  60°.  This  discrepancy  can  perhaps  be 
attributed  to  the  fact  that  an  infinite  strip  is  not  a  very  good  model  for  the  relatively  short 
plate. 
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Figure  7.3:  Triangle-patch  approximation  of  a  vertical,  cylindrical  antenna  that  penetrates 
the  interface  between  two  media. 
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Figure  7.4:  Current  distribution  on  an  open-ended  and  closed  tubular  antenna  of  Fig.  7.3 
partially  buried  in  (a)  dry  earth  and  (b)  salt  water. 
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Figure  7.5:  (aj  Magnitude  and  (b)  phase  of  the  current  J(  on  a  horizontal  tube  partially 
buried  in  a  dielectric  medium.  The  current  is  normalized  to  the  incident  magnetic  field  at 
<f>  =  90°  on  the  surface  of  the  tube. 
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Figure  7.6:  (a)  Magnitude  and  (b)  phase  of  the  current  Je  along  the  center  line  of  a  rectan¬ 
gular  plate  partially  buried  in  dry  earth,  for  the  inclination  angle  a  =  30°.  The  current  is 
normalized  to  the  incident  magnetic  field  at  to  =  0.125  m. 

9 


Normalized  phase  of  Jj  [degree 


Figure  7.9:  Magnitude  of  (a)  the  longitudinal  and  (b)  transverse  component  of  the  current 
induced  on  the  rectangular  plate  of  Fig.  7.6  for  a  =  60°,  by  a  normally-incident  plane  wave 


7.2  Thin  Wires 


The  next  problem  considered  is  that  of  a  straight,  inclined,  thin-wire  antenna  which  is 
partially  buried  in  moist  earth,  as  shown  in  Fig.  7.10.  The  computed  current  distributions 
on  the  antenna  for  a  =  45°  and  a  =  80°  are  compared  in  Figs.  7.11a  and  7.116,  respectively, 
with  the  NEC  [27].  The  agreement  is  seen  to  be  excellent  in  the  first  ca,e,  but  it  is  poor  in 
the  last  case.  A  possible  explanation  of  this  discrepancy  of  the  results  for  a  =  80°  has  been 
given  in  the  previous  section. 

Next,  we  consider  the  case  of  a  vertical,  rectangular,  loop  antenna  partially  immersed  in 
water,  as  illustrated  in  Fig.  7.12,  for  which  measured  data  are  available  [78].  The  current 
distributions  on  the  lower  arm  of  the  loop  are  presented  in  Fig.  7.13  (er  =  81,  a  =  0)  for 
tap  water,  and  in  Figs.  7.14  (er  =  79,  a  —  lS/m)  and  7.15  (er  =  7 6,  cr  =  1.75  S/m)  for 
salt  water.  One  observes  a  good  agreement  of  the  computed  and  measured  results,  both  in 
magnitude  and  phase. 

7.3  Coax— Fed  Microstrip  Patch  Antennas 

In  this  section,  we  present  sample  computed  and  measured  input  impedance  data  (normalized 
to  50  fl)  for  triangular,  rectangular,  and  square  coaxially  fed  microstrip  patch  antennas,  as 
illustrated  in  Fig.  5.7.  In  all  cases,  the  substrate  parameters  are:  er=2.484  and  tan  6  = 
6  x  10-4,  and  the  dimensions  of  the  probe  are:  a  =  0.635  mm  and  6  =  2.095  mm. 

The  results  in  Fig.  7.16  are  for  a  triangular  patch  antenna  previously  analyzed  by  Pi- 
chon  et  al.  [43]  using  a  simple,  zero-order  coax  feed  model.  In  the  numerical  method,  the 
conducting  patch  was  modeled  by  121  triangular  elements.  Computed  results  are  presented 
for  both  the  1st-  and  2nd-order  models  of  the  coax  feed  (cf.  Section  5.3).  As  expected,  the 
2nd-order  model  data  are  closer  to  the  measured  results  than  the  simplified  model  data.  In 
Fig.  7.17  are  shown  results  for  the  same  triangular  patch  antenna,  but  driven  at  a  vertex. 
This  case  illustrates  the  flexibility  of  the  rigorous  coax-feed  model  presented  in  Section  5.3. 
In  Figs.  7.18  and  7.19,  we  present  results  for  a  rectangular  patch  antenna  with  two  different 
substrate  thicknesses.  This  antenna  was  previously  analyzed  by  Hall  and  Mosig  [46],  who 
used  a  rectangular  element  model  of  the  patch.  In  the  present  case,  the  antenna  was  ay  - 
proximated  by  160  triangular  elements.  The  results  indicate  that,  as  expected,  the  lst-order 
feed  model  breaks  down  for  thick  substrates. 

Finally,  in  Fig.  7.20,  we  show  results  for  a  square  patch  antenna  fed  by  a  coaxial  probe 
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Figure  7.11:  Current  distribution  on  the  thin-wire  antenna  of  Fig.  7.10  for  (a)  a  =  45°  and 
(b)  q  =  80°. 
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Figure  7.13:  (a)  Magnitude  and  (b)  phase  of  the  current  J(  on  the  lower  arm  of  the 
rectangular-loop  antenna  of  Fig.  7.12  (er  =  81,  cr  =  0).  The  current  is  normalized  to 
its  value  at  the  point  =  0. 
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Figure  7.14:  (a)  Magnitude  and  (b)  phase  of  the  current  Jt  on  the  lower  arm  of  the 
rectangular-loop  antenna  of  Fig.  7.12  (er  =  79,  a  =  1  S/m).  The  current  is  normalized 
to  its  value  at  the  point  Iq  =  23.3  cm. 


Figure  7.15:  (a)  Magnitude  and  (b)  phase  of  the  current  Jt  on  the  lower  arm  of  the 
rectangular-loop  antenna  of  Fig.  7.12  (er  =  76,  <r  =  1.75  S/m).  The  current  is  normal¬ 
ized  to  its  value  at  the  point  i0  =  23.3  cm. 
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at  an  edge.  The  patch  was  approximated  by  128  triangular  elements.  As  in  the  other  cases' 
presented  above,  we  observe  a  good  agreement  between  the  measured  data  and  the  results 
computed  using  the  2nd-order  coax-feed  model,  except  for  a  small  phase  shift,  which  we 
attribute,  in  part,  to  measurement  inaccuracies  and  to  discrepancies  between  the  catalog 
and  actual  values  of  er. 

7.4  Open  Microstrip  Transmission  Lines 

In  this  section  we  present  results  for  open  microstrip  lines  of  the  form  shown  in  Fig.  5.9, 
except  that  the  conductor  is  confined  to  the  upper  medium.  The  current  distributions 
included  here  are  normalized  to  have  a  maximum  value  of  one  for  the  longitudinal  current 
density. 

The  first  set  of  results  is  for  an  essentially  planar  microstrip  line,  illustrated  in  Fig.  7.21. 
The  PEC  strip  is  of  width  w  and  it  may  be  infinitely  thin,  or  it  may  have  finite  thickness 
t.  In  Fig.  7.22,  we  present  the  dispersion  curves  for  the  lowest  mode  (EHo)  and  the  first 
higher  mode  (EHj)  for  an  infinitely  thin  microstrip  line  and  for  a  line  with  finite  thickness 
( t/w  =  0.1).  The  dimensions  are  w  =  15mm,  d  =  0.794mm,  tT  —  2.32.  In  this  figure,  Oliner 
and  Lee’s  results  [73,  79]  for  t  =  0  are  also  plotted  for  comparison.  Oliner  and  Lee’s  analysis 
is  based  on  the  transverse-resonance  method  in  conjunction  with  the  Wiener-Hopf  approach 
developed  in  [80].  Our  result  for  the  phase  constant  of  the  EHi  mode,  shown  in  Fig.  7.22a,  is 
seen  to  agree  completely  with  Oliner  and  Lee’s  work  (within  the  error  in  reading  from  their 
curves)  both  in  the  bound  and  leaky  regimes,  except  at  the  low  end  of  the  frequency  range. 
The  agreement  of  the  attenuation  constants,  plotted  in  Fig.  7.226,  appears  to  be  slightly  less 
favorable.  Note  that  in  [73]  and  in  [79],  data  for  a  are  only  given  foi,  respectively,  /  >  5  GHz 
and  /  >  6  GHz.  Iri  Figs.  7.18  and  7.19,  we  show  the  current  distributions  of  the  EHo  and 
EHt  modes  on  an  infinitely  thin  strip  at  /  =  5  GHz  and  /  =  10  GHz,  respectively.  The 
longitudinal  currents  (J,)  are  plotted  in  Figs.  7.23a  and  7.24a,  and  the  transverse  currents 
(Jr)  in  Figs.  7.236  and  7.246.  The  symbols  correspond  to  the  locations  where  the  current 
values  are  actually  computed.  Observe  that  nonuniform-width  basis  functions  were  used  to 
better  capture  the  singular  behavior  of  the  longitudinal  current  near  the  edges  of  the  strip. 
It  is  noted  in  Figs.  7.23  and  7.24  that  there  is  very  little  change  in  the  longitudinal  current 
distribution  as  the  mode  passes  from  the  bound  regime  to  the  leaky  regime.  However,  there 
is  a  noticeable  change  in  the  transverse  current  distribution. 

In  Figs.  7.25  and  7.26,  we  show  the  current  distributions  of  the  EHo  and  EHi  modes  for 


Figure  7.21:  Open  microstrip  line. 
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Figure  7.22:  Variation  with  frequency  (a)  of  the  normalized  phase  constant  for  the  lowest 
mode  (EHo)  and  the  first  higher  mode  (EHi),  and  (b)  of  the  normalized  attenuation  constant 
for  the  EH\  mode. 
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Figure  7.23:  (a)  Longitudinal  and  (b)  transverse  current  distributions  of  the  EHo  mode  at 
/  =  5  GHz  and  /  =  10  GHz  for  an  infinitely  thin  microstrip. 
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Figure  7.25:  (a)  Longitudinal  and  (b)  transverse  current  distributions  of  the  EHo  mode  at 
/  =  5  GHz  and  /  =  10  GHz  for  a  microstrip  line  of  finite  thickness. 
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a  strip  with  finite  thickness  (t/w  =  0.1)  at  /  =  5 GHz  and  /  =  10GHz,  respectively.  It  is 
noted  that  the  current  is  concentrated  on  the  lower  face  of  the  strip,  as  expected. 

In  Fig.  7.27,  we  present  the  dispersion  curves  for  the  lowest  mode  (EHo)  and  the  first 
three  higher  modes  (EHi,  EH2,  and  EH3)  for  an  infinitely  thin  microstrip  line  with  w  — 
3  mm,  d  =  0.635mm,  and  tT  —  9.8.  One  should  keep  in  mind  that  the  actual  boundary 
corresponding  to  the  leaky  mode  region  Hi  for  each  mode  is  (1  +  (a/ko)2)* ,  which  obviously 
depends  on  the  attenuation  constant  of  the  mode  (c/.  (6.44)).  As  is  evident  in  Fig.  7.27,  the 
phase  constants  for  the  higher  modes  increase  again  after  reaching  a  minimum  and  continue 
to  increase  as  the  frequency  is  lowered  further.  It  is  noted  that  kv  stays  in  region  Hi  as  the 
frequency  decreases.  In  Fig.  7.28,  we  show  the  normalized  phase  constant  for  the  EH!  mode 
and  the  region  Hi  boundary  corresponding  to  this  mode.  In  Fig.  7.29,  we  present  the  current 
distribution  for  the  EH2  mode  at  /  =  25  GHz  (leaky  regime  in')  and  /  =  35  GHz  (bound 
regime).  In  Fig.  7.30,  we  show  the  current  distribution  for  the  EH3  mode  at  /  =  35 GHz 
(leaky  regime  Hi)  and  at  /  =  50 GHz  (leaky  regime  ii). 

In  Fig.  7.31,  we  plot  the  effective  dielectric  constant,  defined  as  eeg  =  (/ 3/ko )2,  as  a 
function  of  d/ A0,  for  various  widths  and  thicknesses  of  the  strip.  One  observes  that  the 
effective  dielectric  constant  decreases  with  the  thickness  of  the  strip.  This  effect  is  more 
pronounced  at  low  frequencies. 

We  should  mention  that,  at  this  writing,  the  results  for  an  open  microstrip  line  with  finite 
thickness  by  a  rigorous  approach  are  not  available.  We  only  can  compare  our  results  with 
those  for  a  shielded  microstrip  with  finite  thickness.  The  result  in  Fig.  7.31  has  been  found 
to  agree  favorably  with  that  given  in  [81]  (which  is  not  shown  in  the  figure)  for  a  shielded 
strip  line  with  a  large  size  of  the  outer  shield. 

The  last  set  of  results  is  for  a  circular-wire  transmission  line  backed  by  a  grounded 
dielectric  slab,  for* which  a  limited  amount  of  data  are  available  in  the  literature  [82].  In 
Fig.  7.32,  we  present  the  fundamental  mode  effective  dielectric  constant  as  a  function  of  dj -V 
The  corresponding  data  taken  from  [82]  are  also  plotted  for  comparison  and  are  seen  to  agree 
well  with  our  results.  In  Fig.  7.33,  we  present  dispersion  curves  of  the  fundamental  mode 
(EHo)  and  the  first  higher  mod**  (EH\)  for  a  circular-wire  transmission  line  with  h/d= 0.25. 
It  is  observed  that  in  the  chosen  frequency  range  the  EH\  mode  is  leaky  troughout.  Finally, 
in  Fig.  7.34  we  show  the  longitudinal  and  transverse  modal  current  distributions  for  both 
modes  at  <f/Ao=0.3. 
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Figure  7.26:  (a)  Longitudinal  and  (b)  transverse  current  distributions  of  the  EH i  mode  at 
/  =  5  GHz  (leaky  regime  iii)  and  /  =  10  GHz  (bound  regime  i)  for  a  microstrip  transmission 
line  of  finite  thickness. 
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Figure  7.27:  Variation  with  frequency  (a)  of  the  normalized  phase  constant  for  the  lowest 
mode  (EHo)  and  the  first  three  higher  modes  (EHi,EH2,EHs),  and  (b)  of  the  normalized 
attenuation  constants  for  the  EH\,  EH?,  and  EH$  modes  of  am  infinitely  thin  microstrip 
transmission  line. 


Figure  7.29:  (a)  Longitudinal  and  (b)  transverse  current  distributions  of  the  EHi  mode 
at  /  =  25  GHz  (leaky  regime  Hi)  and  /  =  35  GHz  (bound  regime  t)  for  an  infinitely  thin 
microstrip  line. 
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Figure  7.33:  Dispersion  curves  of  the  fundamental  mode  and  the  first  higher  mode  for  a 
circular-wire  transmission  line,  (a)  Phase  constants,  (b)  Attenuation  (leakage)  constant  for 
the  higher  mode. 


Figure  7.34:  (a)  Longitudinal  and  (b)  transverse  current  distributions  of  the  fundamental 
mode  and  the  first  higher  mode  on  a  circular-wire  transmission  line. 


122 


Chapter  8 
Conclusions 


In  this  work,  a  rigorous  and  general  procedure  has  been  developed  for  the  analysis  of  radi¬ 
ation,  scattering,  and  guidance  of  electromagnetic  fields  by  conducting  objects  of  arbitrary 
shape  embedded  in  a  medium  consisting  of  an  arbitrary  number  of  planar,  dielectric  layers. 
The  key  step  in  this  procedure  is  the  transformation  of  the  electric  field  integral  equation  into 
a  mixed-potential  form,  which  is  amenable  to  the  well-established  numerical  solution  tech¬ 
niques  originally  developed  for  arbitrarily-shaped  objects  in  free  space.  Three  particularly 
useful  mixed-potential  integral  equations  (MPIEs)  are  derived  and  their  properties  discussed. 
One  of  the  three  MPIEs,  called  Formulation  C,  which  is  found  to  be  especially  well  suited  for 
the  application  of  the  moment  method,  is  implemented  to  analyze  arbitrarily-shaped,  open 
or  closed,  conducting  surfaces,  which  penetrate  the  interface  between  contiguous  dielectric 
half-spaces.  Thin-wire  structures  are  treated  as  special  cases.  Formulation  C  is  also  special¬ 
ized  to  the  case  of  an  open  transmission  line  consisting  of  an  infinite  conductor  of  arbitrary 
cross-section  partially  embedded  in  a  grounded  dielectric  slab. 

Since  at  this  writing  the  capability  of  analyzing  the  problems  of  electromagnetic  radiation 
and  scattering  by  three-dimensional  PEC  surfaces  of  arbitrary  shape  in  layered  media  does 
not  exist,  the  task  of  validating  our  results  was  a  difficult  one.  Data  for  comparison  were 
only  available  for  thin-wire  antennas,  the  ground  stake  antenna,  and  for  microstrip  patch 
antennas  of  several  simple  shapes.  In  other  cases,  we  had  to  rely  on  the  results  of  two- 
dimensional  analyses.  With  a  few  exceptions,  our  results  presented  in  Chapter  7  compare 
favorably  with  the  few  available  measured  and  numerical  results,  thus  demonstrating — we 
believe — the  validity  of  the  approach  advanced  here. 

Although  open  microstrip  lines  have  been  analyzed  by  both  spectral-domain  [83,  84]  and 
space-domain  [85,  86]  integral  equation  methods,  to  our  knowledge  no  results  have  been 
published  for  microstrips  of  finite  thickness.  Even  for  infinitely  thin  strips,  results  are  only 


available  for  the  bound  modes.  This  is  possibly  due  to  the  fact  that  the  integration  paths 
of  the  spectral-plane  integrals  that  constitute  the  kernels  of  the  integral  equations  must  be 
chosen  differently  for  the  bound  and  leaky  modes,  respectively,  as  explained  in  Section  6.2. 
Hence,  a  computer  code  that  is  successful  for  bound  mode  is  likely  to  fail  to  find  a  leaky 
mode,  unless  special  precautions  are  taken  in  the  program.  This  prompted  some  researchers 
to  voice  skepticism  about  the  very  existence  of  the  leaky  modes  in  a  microstrip.  We  hope  that 
the  results  presented  in  Chapter  7  will  put  to  rest  the  controversy  regarding  the  existence  of 
leaky  microstrip  modes. 

We  wish  to  point  out  that,  at  present,  matrix  fill  time,  not  the  available  computer  mem¬ 
ory,  is  the  overriding  factor,  which  puts  practical  limits  on  the  size  of  objects  that  can  be 
analyzed  by  the  technique  advanced  here.  This  is  due  to  the  fact  that  the  matrix  elements 
comprise  improper,  Sommerfeld  integrals,  which  must  be  repeatedly  evaluated  by  numerical 
quadrature,  as  is  discussed  in  Chapter  6.  Although  it  is  possible,  and  desirable,  to  develop 
analytical  approximations  for  the  Sommerfeld  integrals-and  thus  drastically  reduce  the  com¬ 
putational  expense-we  have  opted  in  this  work  for  rigorous  treatment  of  those  integrals,  to 
avoid  the  uncertainties  associated  with  approximations,  whose  ranges  of  applicability  are 
usually  not  well-defined.  In  the  case  of  microstrip  patch  antennas,  considerable  savings 
in  matrix  fill-time  can  be  realized  by  using  an  interpolation  method  for  the  Sommerfeld 
integrals  (c/.  Section  5.3). 

Although,  for  simplicity,  we  have  concentrated  in  Chapter  5  on  the  analysis  for  objects 
embedded  in  a  two-layer  medium,  we  could  almost  as  easily  treat  the  n-layer  case  (n  >  2), 
for  which  explicit  MPIEs  are  given  in  Chapter  4.  Also,  the  approach  developed  here  can 
be  used  to  analyze  arbitrary  configurations  of  surfaces  and  wires.  One  could,  for  example, 
easily  adapt  the  JUNCTION  code  [7]  for  that  purpose  -  as  is  demonstrated  in  Section  5.3, 
where  the  attachment  mode  used  in  JUNCTION  was  employed  to  solve  the  problem  of  a 
microstrip  patch  antenna  driven  by  a  coaxial  cable.  This  made  it  possible  to  rigorously 
analyze  microstrip  antennas  coax-fed  at  an  edge  or  vertex  of  the  patch.  The  approach 
developed  here  qan  also  be  directly  used  to  calculate  resonant  frequencies  of  microstrip 
antennas  of  arbitrary-possibly  exotic-shapes.  One  could  also  extend  it  to  analyze  infinite 
and  finite  arrays  of  microstrip  patch  antennas  of  arbitrary  shape. 

Future  work  in  this  area  might  involve  dielectric  bodies  of  arbitrary  shape  in  layered 
media.  For  homogeneous  bodies,  one  could  apply  the  surface  integral  formulation  in  con¬ 
junction  with  the  triangle-patch  code  [6]  and  for  inhomogeneous  objects,  the  volume  integral 
approach  in  conjunction  with  the  tetrahedral  element  model  [10].  These  techniques  can  also 
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be  specialized  to  the  analysis  of  dielectric  waveguides  consisting  of  infinite  dielectric  slabs  of 
arbitrary  cross-section  embedded  in  layered  media. 

Finally,  as  was  already  mentioned,  there  is  room  for  improvements  in  the  evaluation  of 
the  Sommerfeld  integrals,  since  most  of  the  computational  effort  is  spent  on  them. 


Appendix  A 

Kernel  Elements  of  Formulation  C 
for  Contiguous  Half-Spaces  and  for  a 
Grounded  Slab 


A.l  Kernel  Elements  for  Contiguous  Half-Spaces 


The  expressions  for  the  elements  of  the  kernels  in  (5.1 )— f 6.3)  can  be  written  down  by  spe¬ 
cializing  the  formulas  given  in  Section  4.3  to  the  two-layer  case  (Fig.  5.1a).  Thus,  when  the 
source  and  observation  points  are  in  the  ith  region  ( i  =  1,2),  one  obtains 
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Ki;x  =  -K'x\  and  K”  =  -K"z 


(A. 5) 


where  Kyz  is  given  by  (2.4)  with  cos  £  replaced  by  sin£.  In  the  above,  we  have  introduced 
the  notation  (cf.  (3.4)) 
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and  r"  —  r'  +  222:'.  The  index  p  in  (A.1)-(A.8)  assumes  the  values  1  or  2,  but  not  equal  to 
i  (z.e.,  p  =  1, 2  and  p  ^  i). 

When  the  source  point  is  in  region  i  and  the  observation  point  in  region  m  /  i,  one 
obtains 
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K™  =  ~^K™  and  K™  =  --K™ 
Hi  Hi 


(A. 14) 


where  K ™'  is  given  by  (A. 11)  with  cos£  replaced  by  sin(.  In  these  equations,  we  have 
introduced  the  notation 
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Several  comments  are  in  order  concerning  the  form  of  the  above  expressions,  which 
resulted-after  a  considerable  amount  of  manipulations-from  the  corresponding  formulas  de¬ 
rived  in  Section  4.3.  The  chosen  forms  have  the  advantage  that  most  of  the  contribution 
to  the  value  of  a  given  kernel  element  comes  from  the  closed  form  term,  thus  deemphasiz¬ 
ing  the  importance  of  the  Sommerfeld  integral.  We  also  note  that  for  the  given  r  and  r', 
only  three  distinct  Sommerfeld  integrals  are  called  for.  Another  advantage  of  this  particular 
formulation  is  that  the  Sommerfeld  integrals  are  well-behaved,  even  when  r  =  r'  on  the 
interface. 


A. 2  Kernel  Elements  for  a  Grounded  Slab 


We  now  specialize  the  formulas  of  Section  4.3  to  the  grounded  slab  case  (see  Fig.  5.16). 
When  the  source  and  observation  points  are  both  in  region  1,  we  obtain 
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When  the  source  point  is  in  region  1  and  the  observation  point  is  in  region  2  (i  =  1,  m  =  2), 
we  obtain 
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When  the  source  point  is  in  region  2  and  the  observation  point  is  in  region  l(i  =  2,m=l), 
we  obtain 
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Finally,  when  the  source  and  observation  points  are  in  region  2  (i  =  2,  m  =  2),  we  obtain 
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where  z<  and  z>  denote,  respectively,  the  lesser  and  the  greater  of  z  and  z' .  K™  and  K™ 
(m  =  1,2,  i  =  1,2)  are  given  by  K™  and  K™,  respectively,  with  cosC  replaced  by  sin(.  In 
the  above  equations,  we  have  introduced  the  notation 
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